anchor.diff:      end do
       rewind(unt)
   
  -    data_type = 'protein'
  -    if (nnuc > nprot) data_type = 'nuc_acid'
  -
  -    call fasta_init(data_type)
  +    if (nnuc > nprot) then
  +       data_type = 'nuc_acid'
  +    else
  +       data_type = 'protein'
  +    end if
   
       string = ""
       k = 1
  
positive.code:   public :: protein_alphabet
    public :: nuc_acid_alphabet
    public :: fasta_read
    character(len=1) :: protein_alphabet(21) = &
         ['A','C','D','E','F','G','H','I','K',&
         'L','M','N','P','Q','R','S','T','V','W','Y','-']
    character(len=1) :: nuc_acid_alphabet(6) = &
         ['A','C','G','T','U','-']
    character(len=21) :: protein_set='ACDEFGHIKLMNPQRSTVWY-'
    character(len=6)  :: nuc_acid_set='ACGTU-'
    integer          :: aamap(239) ! 239 is n. of ascii codes (0 is null)
  
  contains
  
    subroutine fasta_init(data_type)
      character(len=*), intent(in) :: data_type
      ! set the aa -> class mapping
      integer :: k
      
      aamap = 0
      select case(data_type)
      case('protein')
         do k = 1,21
            aamap(iachar(protein_alphabet(k))) = k
         end do
      case('nuc_acid')
         do k = 1,6
            aamap(iachar(nuc_acid_alphabet(k))) = k
         end do
      end select
      
    end subroutine fasta_init
  
    function is_nuc_acid(string) result(res)
      character(len=*), intent(in)  :: string
      logical :: res
  
      res = .false.
      if (verify(trim(string),trim(nuc_acid_set)) == 0) res = .true.
      
    end function is_nuc_acid
    
    function is_protein(string) result(res)
      character(len=*), intent(in)  :: string
      logical :: res
  
      res = .false.
      if (verify(trim(string),trim(protein_set)) == 0) res = .true.
      
    end function is_protein
  
    function aa_to_class(aa) result(c)
      ! takes an a.a. as input, return its index
      character(len=1), intent(in)  :: aa
      integer :: c
      
      c = aamap(iachar(aa))
  
    end function aa_to_class
  
    subroutine fasta_string2seq(string,seq,error_code,error_string)
      ! read a sequence, return an array of integer
  
negative.code:     integer,          intent(out) :: seq(:)
      integer,          intent(out) :: error_code
      character(len=*), intent(out) :: error_string
      integer :: i,n,ind
      integer :: err
  
      error_string = ''
      n = len_trim(string)
      seq = [(index(protein_set,string(i:i)), i=1,len_trim(string))]
      
    end subroutine fasta_string2seq
  
  