c
c                       Los Alamos National Laboratory
c
c       Copyright, 1989.  The Regents of the University of California.
c       This software was produced under a U.S. Government contract
c       (W-7405-ENG-36) by Los Alamos National Laboratory, which is
c       operated by the University of California for the U.S. Department
c       of Energy.  The U.S. Government is licensed to use, reproduce,
c       and distribute this software.  Permission is granted to the
c       public to copy and use this software without charge, provided
c       that this Notice and any statement of authorship are reproduced
c       on all copies.  Neither the Government nor the University makes
c       any warranty, express or implied, or assumes any liability or
c       responsibility for the use of this software.
c
c
c seqf1, CFT77 Fortran version, based on VAXTRAN version
c
c
c CFT77: VAX, SUN, Multiflow dbd
c
c modified (11/08/83) by M.Kanehisa
c modified (12/30/83) by M.Kanehisa
c modified (07/11/84) by C.Burks for longer query sequences.
c     note that, while increasing query length (ndim),
c     the following relations must hold --
c          ldim=2*ndim
c          lwork=.5*ndim
c modified (04/05/87) by C.Burks, added some header comments for clarity
c
c very modified by Dan Davison 10/16/87 for VMS; also handles up to
c  80,000 character sequences.  This should be modified to use the VAX
c  LIB$GET_VM routine (the moral equivalent of the CALL MEMADJ() routine
c  in the CFTLIB library).  My changes are erratically annotated; for best
c  results, use the VMS "diff" program to compare this program with the Cray
c  program of the same name, SEQF1.C (<- this "C" refers to Cray, not the
c  language!)
c
c modified (4/27/88) by D. Davison for Sun, Multiflow
c modified 5/3/89 by dbd for Y/MP & generic ftn 77
c
c   a fortran 77 program to search the dna database for good local
c   homologies.  first the program locates separate diagonal positions
c   containing candidates of good matches by the k-tuple method (wilbur
c   and lipman, proc. natl. acad. sci. usa 80, 726-730, 1973).  then
c   the needleman-wunsch-sellers algorithm (goad and kanehisa, nucleic
c   acids res. 10, 247-263, 1982) is performed to find the best local
c   homology along each of these diagonals.  an option also exists to
c   use the vectorizable nws algorithm for the initial scan of finding
c   the position of the best local homology.
c
c note that all GenBank-like files used as input for this program
c   must be in all-lower-case.
c
c   input parameters
c   ---------------
c   sequence data
c     iseq(1)...iseq(n), sequence to be compared (length n)
c     imin, start value in real numbering scheme
c     iname, sequence identifier
c
c   measure of similarity and output control
c     dmat, weight for a match (default, -1)
c     drep, weight for a mismatch (default, 2)
c     del, weight for a deletion (default, 3)
c     maxd, maximum homology score for each alignment to be printed
c     ktpl, k-tuple size for optional fast scan (default, 4)
c     nw, window size for searching along the diagonal (default, 20)
c     lwid, number of characters per line at the terminal (60 to 120)
c
c   - if ktpl is zero, the nws algorithm is used to scan the sequence
c     for finding the location of the best local homology.
c
c   machine dependence
c   ------------------
c   bit and character manipulations (read instructions in program seqh)
c     lbit, number of bits per word
c
c   storage requirements
c     ndim, maximum length of the query sequence
c     mdim, maximum length of the library sequence
c     ldim, size of work areas
c
c   - the total storage area required is the sum of the following:
c       iseq ... ndim
c       kseq ... ndim
c       jseq ... mdim
c       dist ... max( 2*(ndim+1), ldim )
c       path ... max( 2*ldim, ((2*nwmax)/lpat)*ndim )
c
c   - note the following storage areas are equivalenced.
c       nwint  <- next ->
c       ktuple <- dist ->   <- diag -> <- last ->
c       nwacr  <-- dist -->
c       nwvect <---- zist ---->
c       nwband <-- dist --> <------ path ------->
c
c   - the size of mdim can be dynamically adjusted in cdc and cray.
c
c   files used
c   ----------
c     unit 1, input file of query sequences
c     unit 2, input file of library sequences (database to be searched)
c     unit 6, output file of results (default, user's terminal)
c
c-----------------------------------------------------------------------
c
c   database search
      program seqf(tty,input=tty,output=tty)
c     program seqf(tty,input=tty,output=tty,tape7=tty,tape6=tty)
c      program seqf
      implicit integer(a-y)
      parameter (ndim=2000,mdim=8000,ldim=ndim*2,ktmax=4,nwmax=100)
      character iname*10,jname*70,defin*80,ifile*10,jfile*10,ans*1
      character list(100)*10
      character  kfile*80
      character nucl(5), cmpl(5)
      integer dis(100),loc(100)
      logical logvar
*     level 2,//
      common jseq(mdim)
      common /metric/dmat,drep,del,imin1,lwid
      common /bits/lpat,lres
      common /base/nucl,cmpl
      common /info/jname,defin
      data dmat,drep,del/-1,2,3/,nw/20/
      data nucl/'a','c','g','t','x'/,cmpl/'t','g','c','a','x'/
***
*   vax, multiflow, sun
*      data lbit/32/
*   cdc7600
*     data lbit/60/
*   cray-1
      data lbit/64/
***
      lsize=mdim
      lpat=lbit/3
      lres=lbit-lpat*3
      print *,' '
      print *,'los alamos sequence analysis system -- dna data bank sear
     .ch'
      print *,' '
    1 print 101
  101 format(' input file or help ')
      read 102,ifile
  102 format(a)
      if(ifile.eq.'help') go to 91
c  dd 10/15/87
      inquire(file=ifile,exist=logvar)
      if(.NOT.logvar) then
       print 220, ifile
       goto 92
      else
       open(unit=1,file=ifile,status='UNKNOWN')
      endif
      read(1,102) iname
      if(iname.eq.'locus'.or.iname.eq.'LOCUS')  then
       lf=1
       rewind 1
      else
       lf=2
       rewind 1
      endif
    2 print 103
  103 format(' library file or help ')
      read 102,jfile
      if(jfile.eq.'help') go to 93
c if queryfile ne jfile, open the library file
      if(jfile.eq.ifile) then
        goto 3
c       call addlun(2,1)
       else
        logvar = .TRUE.
        inquire(file=jfile,exist=logvar)
          if(.NOT.logvar) then
            goto 94
          else
            open(unit=2,file=jfile,status='UNKNOWN')
          endif
      endif
    3 print *,' output file ([cr] for your terminal) '
      kfile=''
      read 102,kfile
      if(kfile.ne.'') then
c dd 9/29/87
c       call addlun(6,7)
        open(unit=7,file=kfile, status='new')
c      else
c       call assign(6,kfile,-1)
      endif
      maxd=-15
      ktpl=4
      lwid=70
c added by dan davison, at the request of Christian Burks
      print 106,maxd,ktpl,nw,lwid,dmat,drep,del
c end addition.  just prints out the values of the alignment parameters
c before asking for changes!
  995   print 105
  105 format(' change parameters? (y/n/h) ')
      read 102,ans
      if(ans.eq.'y') go to 4
      if(ans.eq.'n') go to 5
      if(ans.ne.'h') go to 995
      print 106,maxd,ktpl,nw,lwid,dmat,drep,del
  106 format(/' maxd, alignments with more negative scores are printed (
     .current value:',i4,')'/
     .' ktpl, k-tuple size for initial rapid scan (current valuue:',i2,
     .')'/
     .' nw,   window size of diagonal band (current value:',i3,')'/
     .' lwid, terminal line width (current value:',i4,')'/
     .' dmat, negative weight for a match (current value:',i3,')'/
     .' drep, positive weight for a mismatch (current value:',i2,')'/
     .' del,  positive weight for a deletion (current value:',i2,')'/)
      go to 3
    4 print 107
  107 format(' enter maxd,ktpl,nw,lwid,dmat,drep,del ')
      read *,maxd,ktpl,nw,lwid,dmat,drep,del
      if(ktpl.gt.ktmax) ktpl=ktmax
      if(nw.gt.nwmax) nw=nwmax
      lwid=(lwid/10)*10
      if(lwid.lt.60) lwid=60
      if(lwid.gt.120) lwid=120
    5 print 109,maxd,ktpl,nw,lwid,dmat,drep,del
  109 format(/' parameters are:'/3x,'maxd =',i4,
     . 4x,'ktpl =',i2,4x,'nw =',i3,4x,'lwid =',i5/
     . 3x,'dmat =',i3,4x,'drep =',i3,4x,'del =',i3/)
   11 print 111
  111 format(' query sequence ')
      iname='end'
      read 102,iname
      if(iname.eq.'end'.or.iname.eq.'END') goto 99
      call seqin(lf,iname,imin,n,jseq,*11)
      if(n.gt.ndim) go to 95
      n1=n+1
      imin1=imin-1
      imax=imin1+n
      print 110,iname,imin,imax,n
      if(kfile.ne.'') write(7,110) iname,imin,imax,n
      rewind 2
      kount=0
      length=0
      best=0
      call nwint(n,jseq,ktpl)
      if(ktpl.ne.0) go to 31
   21 call genin(lsize,m,jseq,jname,defin,*41)
      kount=kount+1
      length=length+m
***
*   non-vectorized nws algorithm VAX, multiflow
*       call nwacr(m,jseq,dmin,ii,jj)
*   vectorized nws algorithm
      call nwvect(m,jseq,dmin,ii,jj)
***
      if(dmin.lt.best) then
       best=dmin
       l=1
       list(1)=jname(13:22)
      else if(dmin.eq.best.and.l.lt.100) then
       l=l+1
       list(l)=jname(13:22)
      endif
      if(dmin.gt.maxd) go to 21
      call nwband(ii,jj,jseq,dmin,dmin,maxd,nw)
      go to 21
   31 call genin(lsize,m,jseq,jname,defin,*41)
      kount=kount+1
      length=length+m
      call ktuple(m,jseq,dmin,maxd,nt,dis,loc,nw)
      if(dmin.lt.best) then
       best=dmin
       l=1
       list(1)=jname(13:22)
      else if(dmin.eq.best.and.l.lt.100) then
       l=l+1
       list(l)=jname(13:22)
      endif
      if(nt.le.0) go to 31
      do 32 i=1,nt
       ij=loc(i)
       if(ij.le.m) then
        ii=n
        jj=ij
       else
        ii=n+m-ij
        jj=m
       endif
   32  call nwband(ii,jj,jseq,-999999,dis(i),maxd,nw)
      go to 31
   41 print 120,length,kount
      if(ktpl.eq.0) then
       print 130,best,(list(i),i=1,l)
       if(kfile.ne.'') write(7,130) best,(list(i),i=1,l)
      else
       print 140,best,(list(i),i=1,l)
       if(kfile.ne.'') write(7,140) best,(list(i),i=1,l)
      endif
   43 print 112
  112 format(/' continue? (y/n)')
      read 102,ans
      if(ans.eq.'y') go to 11
      if(ans.eq.'n') go to 99
      go to 43
   91 print 200
      go to 1
   92 print 220,ifile
      go to 1
   93 print 210
      go to 2
   94 print 230,jfile
      go to 2
   95 print *,'*** query sequence too long.'
      go to 11
   99 stop
  110 format(' query sequence',a12,'  from',i8,'  to',i8,'  total',i7)
  120 format(/' total of ',i8,' bases in ',i5,' sequences searched')
  130 format(/' the best homology score ',i6,
     .'  was found in:'/5(a13))
  140 format(/' the best homology score (uncorrected) ',i6,
     .'  was found in:'/5(a13))
  200 format(/' input file of query sequences must be either one of the
     .following:'/
     .'   1. genbank format (if the first line is a locus line)'/
     .'   2. stanford seq format (otherwise)'/)
  210 format(/' choose one of the following library files:'/
     .'   dna   entire genbank database'/
     .'   dnam  mammalian sequences'/
     .'   dnan  non-mammalian vertebrate sequences'/
     .'   dnai  invertebrate sequences'/
     .'   dnap  plant sequences'/
     .'   dnao  organelle sequences'/
     .'   dnab  bacteria sequences'/
     .'   dnas  structural rna sequences'/
     .'   dnav  virus sequences'/
     .'   dnaf  phage sequences'/
     .'   dnar  recombined or synthetic sequences'/
     .' or any file in the genbank format.'/)
  220 format(' *** file',a12,' does not exist. try again.')
  230 format(' *** file',a12,' does not exist as library. try again.')
      end
c   initialization for nws algorithm
      subroutine nwint(m,jseq,ktpl)
      implicit integer(a-y)
      parameter (ndim=2000,mdim=8000,ldim=ndim*2,ktmax=4,nwmax=100)
      parameter (lwork=250)
*     level 2,jseq
      character*1  blank, hyphen, colon
      character  kfile*80
      character nucl(5), cmpl(5)
      dimension iseq(ndim),jseq(mdim),dist(ndim*2+2),diff(5,5)
      dimension zist(ndim*3+3)
      dimension kseq(ndim),next(ndim),mark(4**ktmax+1)
      dimension diag(ldim),last(ldim),work(2,lwork),dis(1),loc(1)
c      dimension pat(nwmax*2),path(ldim*2)
      integer pat(nwmax*2),path(ldim*2)
      dimension ishf(13),jshf(13),inum(13),jnum(13),midl(120)
      equivalence (dist(1),next(1),zist(1))
      equivalence (path(1),diag(1)),(path(ldim+1),last(1))
      common /metric/ddmat,ddrep,ddel,imin1,lwid
      common /bits/lpat,lres
      common /base/nucl,cmpl
      common /info/jname,defin
      character jname*70,defin*80,fm*17,lchr*2
      logical final,first
      save n,n1,iseq,diff,dmat,drep,del,dmin,dmkt
      save kseq,len,kt,kt1,first
      data blank,hyphen,colon/' ','-',':'/
      iand(i,j)=i.and.j
      ior(i,j)=i.or.j
c      shift(i,j)= ishft(i,j)
      ishft(i,j)=shift(i,j)
c      print *,' statement functions made'
      lpat3=lpat*3
      lc=lwid-1
      dmat=ddmat
      drep=ddrep
      del=ddel
      zmat=dmat
      zrep=drep
      zel=del
c      print *,' constants set in nwint'
      do 12 i=1,4
        do 10 j=1,4
   10     diff(i,j)=drep
        diff(i,5)=0
        diff(5,i)=0
   12   diff(i,i)=dmat
      diff(5,5)=0
      n=m
      do 14 i=1,n
   14   iseq(i)=jseq(i)
      n1=n+1
      if(ktpl.eq.0) return
c   initialization for k-tuple search
c      print *,' initializing k-tuple search'
      do 20 i=1,n
   20   kseq(i)=0
      kt=ktpl
      kt1=kt-1
      dmkt=dmat*kt
      len=4**kt
      do 22 k=1,len+1
   22   mark(k)=0
      do 28 i=1,n-kt1
        i1=i-1
        k=0
        do 24 j=1,kt
          if(iseq(i1+j).eq.5) go to 25
   24     k=k*4+iseq(i1+j)-1
        go to 26
   25   k=len
   26   k=k+1
        if(mark(k).eq.0) then
          mark(k)=i
        else
          kseq(next(k))=i
        endif
   28   next(k)=i
c      print *,' exiting ktuple init'
      return
c   k-tuple search
      entry ktuple(m,jseq,ddmin,maxd,nt,dis,loc,nw)
c      print *,' entry into ktuple'
      first=.true.
      msafe=-6
      limit=ldim-n+1
      lwork1=lwork/(m/limit+1)
      nm1=n+m-1
      do 30 i=1,min0(ldim,nm1)
        dist(i)=0
        diag(i)=0
   30   last(i)=0
      mk=m-kt1
      mt=0
      ja=1
   31 ja1=ja-1
      jb=ja+limit-1
      if(jb.gt.mk) then
        jb=mk
        limit=nm1-ja1
      endif
c      print *,' beginning loop 34 in nwint'
      do 34 j=ja,jb
        j1=j-1
        jkt=j-kt
        k=0
        do 32 i=1,kt
   32     k=k*4+jseq(j1+i)-1
        if(k.gt.len) go to 34
        k=k+1
        if(mark(k).eq.0) go to 34
        i=mark(k)
   33   ij=j-i+n-ja1
        if(last(ij).eq.0) then
          dist(ij)=dmkt
c      print *,' calc distance A'
        else
          if(jkt-last(ij).le.0) then
            dist(ij)=dist(ij)+dmat
c      print *,' calc distance b'
          else
            dist(ij)=min0(dist(ij)+(jkt-last(ij))*drep,0)+dmkt
c      print *,' calc distance c'
          endif
        endif
        if(dist(ij).lt.diag(ij)) diag(ij)=dist(ij)
        last(ij)=j
        if(kseq(i).eq.0) go to 34
        i=kseq(i)
        go to 33
   34 continue
   35 l=0
c      print *,' starting limit check'
      do 36 i=1,limit
        if(diag(i).lt.msafe) then
          l=l+1
          if(l.gt.lwork1) go to 39
          work(1,l+mt)=diag(i)
          work(2,l+mt)=i+ja1
        endif
   36 continue
      mt=mt+l
      if(jb.ge.mk) go to 41
      ll=ldim-limit
      do 37 i=1,ll
        dist(i)=dist(i+limit)
        diag(i)=diag(i+limit)
   37   last(i)=last(i+limit)
      do 38 i=ll+1,ldim
        dist(i)=0
        diag(i)=0
   38   last(i)=0
      ja=ja+limit
      go to 31
   39 msafe=msafe-1
*     print *,'*** work area overflow. recovering...'
      go to 35
c   locate good alignments on separate diagonals
c      print *,' checking diagonals'
   41 if(mt-1) 42,43,44
   42 nt=0
      ddmin=0
      return
   43 nt=1
      dis(1)=work(1,1)
      loc(1)=work(2,1)
      ddmin=dis(1)
      return
   44 nt=0
      nw2=nw/2
      ddmin=0
      do 48 j=1,100
      k=0
      dmin=0
      do 45 i=1,mt
        if(work(1,i).ne.0.and.work(1,i).lt.dmin) then
c      print *,' checking work arrays a'
          dmin=work(1,i)
          k=i
        endif
   45 continue
      if(k.eq.0) return
      work(1,k)=0
      ij=work(2,k)
      do 46 i=1,mt
        if(work(1,i).ne.0) then
          if(iabs(work(2,i)-ij).lt.nw2) then
c      print *,' checking work arrays b'
            dmin=dmin+min0(work(1,i)+iabs(work(2,i)-ij)*del,0)
            work(1,i)=0
          endif
        endif
   46 continue
      if(dmin.lt.ddmin) ddmin=dmin
      if(dmin.le.maxd) then
        nt=nt+1
        dis(nt)=dmin
        loc(nt)=ij
      endif
   48 continue
      return
c   row-wise nws algorithm (non-vectorizable)
      entry nwacr(m,jseq,ddmin,ii2,jj2)
      do 50 i=1,n1+1
   50   dist(i)=0
      jn1=n1
      jn2=0
      dmin=0
      do 54 j=1,m
        js=jseq(j)
        do 52 i=1,n
          dist(i+1+jn1)=min0(dist(i+jn2)+diff(iseq(i),js),
     .                  dist(i+1+jn2)+del,dist(i+jn1)+del,0)
          if(dist(i+1+jn1).lt.dmin) then
            dmin=dist(i+1+jn1)
            i2=i
            j2=j
          end if
   52   continue
        jj=jn1
        jn1=jn2
        jn2=jj
   54 continue
      ddmin=dmin
      ii2=i2
      jj2=j2
      first=.true.
      return
c   anti-diagonal nws algorithm (vectorizable)
      entry nwvect(m,jseq,ddmin,ii2,jj2)
c      print *,' entry into nwvec'
      n2=n1+n1
      n3=n2+n1
      m1=m-1
      nm1=n+m1
      do 60 i=1,n3
   60   zist(i)=0.
      jn=0
      jn1=n2
      jn2=n1
      zmin=0
      do 66 k=1,nm1
        k1=k+1
        ia=max0(1,k-m1)
        ib=min0(k,n)
c        print *, ' only important vector loop in the pgm'
cdir$   ivdep
        do 62 i=ia,ib
          z=cvmgt(zmat,zrep,iseq(i).eq.jseq(k1-i))
          zist(i+1+jn)=amin1(zist(i+jn2)+z,
     .                 zist(i+1+jn1)+zel,zist(i+jn1)+zel,0.)
   62   continue
        i=ismin(ib-ia+1,zist(ia+1+jn),1)+ia
        if(zist(i+jn).lt.zmin) then
          zmin=zist(i+jn)
          i2=i
          j2=k1
        end if
        jn=mod(jn+n1,n3)
        jn1=mod(jn1+n1,n3)
        jn2=mod(jn2+n1,n3)
   66 continue
      i2=i2-1
      j2=j2-i2
      ddmin=zmin
      ii2=i2
      jj2=j2
      first=.true.
      return
c   nws algorithm along the diagonal band
      entry nwband(nn,mm,jseq,ddmin,dd,maxd,nw)
      final=.false.
c        print *, ' entry into nwband'
      nn1=nn+1
      mm1=mm+1
      nw1=nw-1
      nw2=nw+nw
      dist(1)=0
      dist(1+nn1)=0
      do 70 i=1,nw
       dist(i+1)=0
   70  pat(i)=0
      do 72 i=nw+2,nn1
       dist(i)=1000
   72  dist(i+nn1)=1000
      k=1
      jn=-nw
      jn1=nn1
      jn2=0
      dmin=0
      do 78 j=1,mm
       if(j.eq.nn) final=.true.
       js=jseq(mm1-j)
       ia=max0(j-nw,1)
       if(ia.gt.1) dist(ia+jn1)=1000
       ib=min0(j+nw1,nn)
       do 74 i=ia,ib
c       print *, ' distance check in nwband'
        d1=dist(i+1+jn2)+del
        d2=dist(i+jn2)+diff(iseq(nn1-i),js)
        d4=dist(i+jn1)+del
        d=min0(d1,d2,d4,0)
        dist(i+1+jn1)=d
        pat(i-jn)=0
        if(d.eq.0) go to 74
        if(d1.eq.d) pat(i-jn)=ior(pat(i-jn),1)
        if(d2.eq.d) pat(i-jn)=ior(pat(i-jn),2)
        if(d4.eq.d) pat(i-jn)=ior(pat(i-jn),4)
        if(d.lt.dmin) then
         dmin=d
         i2=i
         j2=j
         if(dmin.le.ddmin) final=.true.
        end if
   74  continue
       i=0
   75  patwd=0
       do 76 l=1,lpat
        i=i+1
        if (pat(i).lt.0.or.pat(i).gt.7) go to 76
c        patwd=(patwd.or.ishft(pat(i),lpat3-l*3))
c        print *,' piddling on patwd'
c        patwd=ior(patwd,ishift(pat(i),lpat3-l*3))
        patwd=ior(patwd,ishft(pat(i),lpat3-l*3))
   76   if(i.eq.nw2) go to 77
   77  path(k)=patwd
       k=k+1
       if(i.ne.nw2) go to 75
       if(final) go to 81
       jn=jn+1
       jj=jn1
       jn1=jn2
       jn2=jj
   78 continue
c   trace back
   81 if(dmin.gt.maxd) return
c      print *,' beeginning traceback'
      if(first) then
        write(*,110) jname,defin
        if (kfile.ne.'') write(7,110) jname,defin
        first=.false.
        mt=0
      endif
      imin=nn1-i2
      jmin=mm1-j2
      do 82 i=1,mt
   82   if(work(1,i).eq.imin.and.work(2,i).eq.jmin) return
      mt=mt+1
      work(1,mt)=imin
      work(2,mt)=jmin
      l=0
      i=i2
      j=j2
      kn=(nw2-1)/lpat+1
      k=(j2-1)*kn+(i2-(j2-1-nw)-1)/lpat+1
      ks=(mod(i2-(j2-1-nw)-1,lpat)+1)*3
   83  if(i.eq.0.or.j.eq.0) go to 87
***
*   rightward (negative) shift VAX
*        jump=iand(ishft(path(k),ks-lpat3),7)
*       jump=and(lshift(path(k),ks-lpat3),7)
*   leftward circular shift (cray)
       jump=iand(ishft(path(k),ks+lres),7)
***
       if(iand(jump,2).eq.0) go to 84
c       if(and(jump,2).eq.0) go to 84
       l=l+1
       dist(l)=ichar(nucl(iseq(nn1-i)))
       dist(l+ndim)=ichar(nucl(jseq(mm1-j)))
       i=i-1
       j=j-1
       k=k-kn
       go to 83
   84  if(iand(jump,1).eq.0) go to 85
c   84  if(and(jump,1).eq.0) go to 85
       l=l+1
       dist(l)=ichar(hyphen)
       dist(l+ndim)=ichar(nucl(jseq(mm1-j)))
       j=j-1
       k=k-kn
       ks=ks+3
       if(ks.le.lpat3) go to 83
       ks=3
       k=k+1
       go to 83
   85  if(iand(jump,4).eq.0) go to 87
c   85  if(and(jump,4).eq.0) go to 87
       l=l+1
       dist(l)=ichar(nucl(iseq(nn1-i)))
       dist(l+ndim)=ichar(hyphen)
       i=i-1
       ks=ks-3
       if(ks.ge.3) go to 83
       ks=lpat3
       k=k-1
       go to 83
   87 imin=nn1-i2+imin1
      jmin=mm1-j2
      imax=nn1-i+imin1
      jmax=mm1-j
c   alignment output
      nl=l
      match=0
      mismat=0
      idel=0
      jdel=0
      if(imin.lt.0) then
       inum(1)=(imin/10)*10
       ishf(1)=inum(1)-imin+1
       if(iabs(inum(1)).ge.10**(min0(5,ishf(1))-1)) then
        inum(1)=inum(1)+10
        ishf(1)=ishf(1)+10
       endif
      else
       inum(1)=(imin/10+1)*10
       ishf(1)=inum(1)-imin+1
       if(inum(1).ge.10**min0(5,ishf(1))) then
        inum(1)=inum(1)+10
        ishf(1)=ishf(1)+10
       endif
      endif
      jnum(1)=(jmin/10+1)*10
      jshf(1)=jnum(1)-jmin+1
      if(jnum(1).ge.10**min0(5,jshf(1))) then
       jnum(1)=jnum(1)+10
       jshf(1)=jshf(1)+10
      endif
      i=imin
      j=jmin
      l1=1
   91  l2=l1+lc
       if(l2.gt.nl) l2=nl
       k=1
       ki=1
       kj=1
       do 92 l=l1,l2
        if(dist(l).eq.dist(l+ndim)) then
         midl(k)=ichar(colon)
         match=match+1
        else
         midl(k)=ichar(blank)
        endif
        k=k+1
c dd 4/28/88 for f77
        if(char(dist(l)).eq.hyphen) then
         ishf(ki)=ishf(ki)+1
         idel=idel+1
        else
         i=i+1
         if(i.gt.inum(ki)) then
          ki=ki+1
          inum(ki)=inum(ki-1)+10
          ishf(ki)=10
         endif
        endif
c dd 4/28/88 for f77
        if(char(dist(l+ndim)).eq.hyphen) then
         jshf(kj)=jshf(kj)+1
         jdel=jdel+1
        else
         j=j+1
         if(j.gt.jnum(kj)) then
          kj=kj+1
          jnum(kj)=jnum(kj-1)+10
          jshf(kj)=10
         endif
        endif
   92   continue
       k1=k-1
       ki1=max0(1,ki-1)
       kj1=max0(1,kj-1)
       call numbrp(ishf,inum,ki1)
       write(*,100) (char(dist(l)),l=l1,l2)
       if (kfile.ne.'') write(7,100) (char(dist(l)),l=l1,l2)
       write(*,100) (char(midl(k)),k=1,k1)
       if (kfile.ne.'') write(7,100) (char(midl(k)),k=1,k1)
       write(*,100) (char(dist(l+ndim)),l=l1,l2)
       if (kfile.ne.'') write(7,100) (char(dist(l+ndim)),l=l1,l2)
       call numbrp(jshf,jnum,kj1)
       if(l2.eq.nl) go to 93
       l1=l2+1
       ishf(1)=inum(ki)-i+1
       jshf(1)=jnum(kj)-j+1
       inum(1)=inum(ki)
       jnum(1)=jnum(kj)
       if(imin.lt.0) then
        if(iabs(inum(1)).ge.10**(min0(5,ishf(1))-1)) then
         inum(1)=inum(1)+10
         ishf(1)=ishf(1)+10
        endif
       else
        if(inum(1).ge.10**min0(5,ishf(1))) then
         inum(1)=inum(1)+10
         ishf(1)=ishf(1)+10
        endif
       endif
       if(jnum(1).ge.10**min0(5,jshf(1))) then
        jnum(1)=jnum(1)+10
        jshf(1)=jshf(1)+10
       endif
       go to 91
   93 mismat=nl-(match+idel+jdel)
      z=real(match)/real(nl)*100.
      write(*,120) dmin,dd,z,match,nl
       if (kfile.ne.' ') write(7,120) dmin,dd,z,match,nl
      return
  100 format(120a1)
  110 format(/a70/a80/)
  120 format(/6x,'homology score =',i6,4x,'(',i5,' )'/
     . 6x,'percent match  =',f6.1,'%   (',i4,' /',i4,' )'/)
      end
c   variable format
      subroutine numbrp(ishf,inum,ki)
      integer ishf(1),inum(1)
      character line(1)
      character fm*52,lchr*2
      fm(1:4)='('
      l=2
      do 10 k=1,ki
       write(unit=lchr,fmt='(i2)') ishf(k)
       fm(l:l+3)='i'//lchr//','
   10  l=l+4
      fm(l-1:l-1)=')'
      write(*,fm) (inum(k),k=1,ki)
      return
      end
c   genbank input
      subroutine genin(lsize,n,seq,locus,defin,*)
      parameter (ndim=2000,mdim=8000,ldim=ndim*2,ktmax=4,nwmax=100)
      character text*80,locus*70,defin*80, standrd*80
*     level 2,seq
      integer seq(mdim)
      character  kfile*80
      character nucl(5), cmpl(5)
      character alphseq(ndim)
      common /base/nucl, cmpl
      do 10 i=1,1000
       read(2,100,end=21) text
   10  if(text(1:10).eq.'locus'.or.text(1:10).eq.'LOCUS') goto 11
      go to 21
   11 locus=text
      if(locus(13:22).eq.'end'.or.locus(13:22).eq.'END') goto 21
      read(unit=text,fmt='(22x,i7)') n
      if(n.gt.lsize) then
c       call memadj(n-lsize,ierr)
c       lsize=n
        write(*,120) locus, n, lsize
        n=lsize
      else if(lsize.gt.mdim) then
c       call memadj(mdim-lsize,ierr)
       lsize=mdim
      endif
      read(2,100) standrd
      read(2,100) defin
      do 12 i=1,1000
       read(2,100,end=21) text
   12  if(text(1:10).eq.'origin'.or.text(1:10).eq.'ORIGIN') goto 13
      go to 21
   13 read(2,110) (alphseq(i),i=1,n)
      do 16 i=1,n
c       print *, ' letter read in: ', alphseq(i)
       if(alphseq(i).eq.'u') alphseq(i)='t'
       do 14 j=1,4
   14   if(alphseq(i).eq.nucl(j)) go to 16
       j=5
   16  seq(i)=j
      return
   21 return 1
  100 format(a)
  110 format((9x,6(1x,10a1)))
 120  format(//,4x,' Truncated ',a22,', length ', i6,' avail space= ',
     .     i6,'.')
      end
c   sequence input
      subroutine seqin(lf,locus,min,n,seq,*)
      parameter (ndim=2000, mdim=8000)
      implicit integer(a-z)
*     level 2,seq
      dimension seq(mdim)
c     line(80)
      character*10 locus,name,head,ans*1
      character    nucl(5), cmpl(5)
      character alphseq(ndim), line(80)
      common /base/ nucl, cmpl
      rewind 1
      if(lf.eq.1) go to 13
   11 read(1,100,end=51) name
      if(name.ne.locus) go to 11
      go to 17
   13 read(1,100,end=51) head,name
      if(head.eq.'locus') then
        if (name.eq.locus) then
          goto 15
        else
          goto 13
        endif
      else if(head.eq.'LOCUS') then
        if (name.eq.locus) then
          goto 15
        else
          goto 13
        endif
      else
        goto 13
      endif
c     if(name.ne.locus) go to 13
   15 read(1,100) head
      if(head.eq.'origin') then
        goto 17
      elseif (head.eq.'ORIGIN') then
        goto 17
      else
        goto 15
      endif
   17 print 101
  101 format(' start,end (1,0 for the entire sequence) ')
      min=1
      max=0
      read *,min,max
      min1=min-1
      i=0
      if(lf.eq.1) go to 23
   21 read(1,110,end=25) (line(j),j=1,80)
      do 22 j=1,80
       if(line(j).eq.'1') go to 25
       if(line(j).eq.' ') go to 22
       i=i+1
       if(i.lt.min) go to 22
       alphseq(i-min1)=line(j)
       if(i.eq.max) go to 25
   22  continue
      go to 21
   23 read(1,120) head,(line(j),j=1,60)
      if(head.ne.' ') go to 25
      do 24 j=1,60
       if(line(j).eq.' ') go to 25
       i=i+1
       if(i.lt.min) go to 24
       alphseq(i-min1)=line(j)
       if(i.eq.max) go to 25
   24  continue
      go to 23
   25 max=i
      n=max-min+1
      do 32 i=1,n
       if(alphseq(i).eq.'u') alphseq(i)='t'
       do 30 j=1,4
   30   if(alphseq(i).eq.nucl(j)) go to 32
       j=5
   32  seq(i)=j
   33 print 102
  102 format(' complement? (y/n) ')
      read 110,ans
      if(ans.eq.'n') return
      if(ans.ne.'y') go to 33
   41 min=-max
      do 42 i=1,n
   42  seq(i)=ichar(cmpl(seq(i)))
      n1=n+1
      n2=n/2
      do 44 i=1,n2
       s=seq(i)
       seq(i)=seq(n1-i)
   44  seq(n1-i)=s
c      print *,seq
      return
   51 print 200,locus
      return 1
  100 format(a10,2x,a10)
  110 format(80a1)
  120 format(a1,8x,6(1x,10a1))
  200 format(' *** ',a10,' was not found *** try again ***')
      end