Program PSSMsearch
      IMPLICIT NONE
C***********************************************************************
C  declarations of logical unit numbers or input and output
C***********************************************************************
      INTEGER  STDIN, STDOUT, PRUNIT, SEQUNIT

      PARAMETER  (  STDIN = 5  , STDOUT = 6  )
      PARAMETER  (  PRUNIT = 15  , SEQUNIT = 16  )

C***********************************************************************
C   constants defining file types
C***********************************************************************
      INTEGER      GENBNK, EMBL, NBRF, MOLGN, SIMPL
      PARAMETER ( SIMPL = 0, GENBNK = 1, EMBL = 2, NBRF = 3, MOLGN=4) 
      INTEGER      PROFIL, GCG, FASTA
      PARAMETER (PROFIL=5, GCG=6, FASTA=7)
      INTEGER   NA, PROT
      PARAMETER   (NA=1, PROT=2)
!***********************************************************************
!   Dimensions of profile objects.
!   (If POSITION_MAX isn't big enough, the routine PROFILE_SCORE_READ 
!   will tell you what to increase it to, just before it stops).
!***********************************************************************
      integer acid_num
      integer position_max
      parameter ( acid_num = 23 )
      parameter ( position_max = 200 )
C***********************************************************************
C  Sequence Definitions
C***********************************************************************
      INTEGER SEQDP,FRAMEDP,FOREVDP
      PARAMETER (SEQDP=10000)
      PARAMETER (FRAMEDP=6)
      PARAMETER (FOREVDP=4) 
      INTEGER SEQ(FRAMEDP,SEQDP)
      CHARACTER*1 TSEQ(SEQDP)
      CHARACTER*1 SEQUENCE(SEQDP,FOREVDP)
      CHARACTER*32 LOCUS
      CHARACTER*32 ACCESSION
      CHARACTER*256 DEFN
      INTEGER LENGTH
      INTEGER ALPHA(0:127)
      INTEGER TABDIM
      PARAMETER (TABDIM=32)
      CHARACTER*1 TABLE(0:TABDIM)
C***********************************************************************
C   Profile objects.
C***********************************************************************
      character acid_code(acid_num)
      character conserved(position_max)
      integer entropy(position_max)
      integer gap_extend_percent(0:position_max)
      integer gap_open_percent(0:position_max)
      integer position_num
      integer score(position_max,acid_num)
      integer score2(acid_num)
C***********************************************************************
C    Misc. Variables
C***********************************************************************
      INTEGER SEQFMT,ERROR, guessnap
      integer iberr, i, j

C***********************************************************************
C   OPEN UP SEQUENCE FILE
C***********************************************************************
C      print *, 'Opening up file'
C      call flush(6)      
      OPEN( UNIT = SEQUNIT, FILE = 'MY.SEQ', IOSTAT = IBERR, 
     +      STATUS = 'OLD', ACTION='READ')

      IF( IBERR.NE.0) THEN
         WRITE(STDOUT,*) 'PROBLEM OPENING SEQUENCE FILE'
         STOP 'NUM 1'
      ENDIF
C********************************************************************
C   Determine the file format of the sequence; 
C********************************************************************
C      print *, 'Calling whatfm'
C      call flush(6)
      Call WHATFM(sequnit,stdout,seqfmt,ERROR,.false.,guessnap)
C      print *, 'Fileformat=',seqfmt,' GUESSNAP=',guessnap
C      call flush(6)
      IF(GUESSNAP.EQ.NA) THEN
         call nuclab(ALPHA,TABLE,TABDIM)
      ELSEIF (GUESSNAP.EQ.PROT) THEN
         call protab(ALPHA,TABLE,TABDIM)
      ELSE
         print *, 'INVALID ALPHABET'
         call flush(6)
         stop
      ENDIF   
C      print *, 'Alphabet selected'
C      call flush(6)
C********************************************************************
C     Read in the sequence; Only permit single sequences at this point 
C     in time.
C********************************************************************
      CALL READSEQ(ERROR,SEQUNIT,SEQFMT,SEQ(1,1),
     +        SEQDP,SEQUENCE(1,1),SEQDP,LOCUS,ACCESSION,DEFN,
     +        ALPHA,LENGTH)
C      print *, 'Done calling readseq Length=',length
C      DO 30 i=1,LENGTH
C            write (6,fmt='(A1,$)') SEQUENCE(I,1)
C            if (MOD(I,70).EQ.0) print *, ' '            
C 30   CONTINUE

      IF (GUESSNAP.EQ.NA) THEN
          print *, 'CALLING TRANSLATE'
          CALL TRANSLATE(ERROR,SEQUENCE,FOREVDP,SEQDP,TABLE,TABDIM,
     +        ALPHA,TABDIM,LENGTH)
          CALL DUMP_TRANSLATE(ERROR,SEQUENCE,FOREVDP,SEQDP,LENGTH)
           
      ELSE
           
      ENDIF

C      DO 200 j=1,4
C         PRINT *, 'SEQUENCE[',j,'] Length=',length
C         DO 300 i=1,LENGTH
C            write (6,fmt='(A1,$)') SEQUENCE(I,J)
C            if (MOD(I,70).EQ.0) print *, ' '
C 300     CONTINUE
C 200  CONTINUE

C EOF read in from library file.
      IF (ERROR.EQ.3)  THEN
      ENDIF

!********************************************************************
!  Read in the profile data.
!********************************************************************

      open ( unit = PRUNIT, file = 'profile.txt', status = 'old', 
     &  iostat = IBERR, action = 'READ' )

      if ( IBERR /= 0 ) then
        write ( *, * ) ' '
        write ( *, * ) 'PSSMsearch - Fatal error!'
        write ( *, * ) '  Could not open the profile file.'
        stop
      end if

      call profile_score_read ( acid_code, acid_num, conserved, 
     &  entropy, gap_open_percent, gap_extend_percent, PRUNIT, 
     &  position_max, position_num, score, score2 )

      close ( unit = PRUNIT )

!********************************************************************
!  Print the profile data.
!********************************************************************

      call profile_score_print ( acid_code, acid_num, conserved, 
     &  entropy, gap_open_percent, gap_extend_percent, position_max, 
     &  position_num, score, score2 )

!********************************************************************
!  Compute the local optimal alignment.
!********************************************************************

      call local_align ( length, sequence(1,1) )

      STOP
      END