program main

!*****************************************************************************80
!
!! ps_gg_align_test() tests ps_gg_align().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 December 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'ps_gg_align_test():'
  write ( *, '(a)' ) '  FORTRAN90 version.'
  write ( *, '(a)' ) '  Test PS_GG_ALIGN().'

  call test01 ( )
  call test02 ( )
  call test03 ( )
  call test04 ( )
  call test05 ( )
  call test06 ( )
  call test07 ( )
  call test08 ( )
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'PS_GG_ALIGN_TEST'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop 0
end
subroutine test01 ( )

!*****************************************************************************80
!
!! TEST01 tests PROFILE_SCORE_READ, PROFILE_SCORE_PRINT.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 December 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: acid_num = 23
  integer, parameter :: iunit = 1
  integer, parameter :: position_max = 200

  character acid_code(acid_num)
  character conserved(position_max)
  integer entropy(position_max)
  character ( len = 80 ) file_name
  integer gap_extend_percent(0:position_max)
  integer gap_open_percent(0:position_max)
  integer ios
  integer position_num
  integer score(position_max,acid_num)
  integer score2(acid_num)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  PROFILE_SCORE_READ reads profile scoring '
  write ( *, '(a)' ) '    information.'
  write ( *, '(a)' ) '  PROFILE_SCORE_PRINT prints the information.'

  file_name = 'profile.txt'
!
!  Read the data from the file.
!
  open ( unit = iunit, file = file_name, status = 'old', iostat = ios )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST01 - Fatal error!'
    write ( *, '(a)' ) '  Could not open the data file.'
    return
  end if

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

  close ( unit = iunit )
!
!  Print the data.
!
  call profile_score_print ( acid_code, acid_num, conserved, entropy, &
    gap_open_percent, gap_extend_percent, position_max, position_num, &
    score, score2 )

  return
end
subroutine test02 ( )

!*****************************************************************************80
!
!! TEST02 tests PS_GG_FPQ, PS_GG_FSQ.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 December 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: m = 4
  integer, parameter :: n = 5
  integer, parameter :: lds = m

  character b(n)
  real base
  real e(0:lds,0:n)
  real f(0:lds,0:n)
  real, parameter :: gap_extend = -0.5E+00
  integer gap_extend_percent(0:m)
  real, parameter :: gap_open = -2.0E+00
  integer gap_open_percent(0:m)
  real ge(0:m)
  real go(0:m)
  integer i
  integer j
  integer m1
  integer m2
  integer n1
  integer n2
  integer npath
  integer pathi(m+n+1)
  integer pathj(m+n+1)
  real, external :: ps_score_simple
  real s(0:lds,0:n)
  integer t(0:lds,0:n)

  b(1) = 'A'
  b(2) = 'C'
  b(3) = 'G'
  b(4) = 'C'
  b(5) = 'T'

  gap_open_percent(0) =  95
  gap_open_percent(1) =  90
  gap_open_percent(2) =  75
  gap_open_percent(3) =  50
  gap_open_percent(4) =  25

  gap_extend_percent(0) = 100
  gap_extend_percent(1) =  80
  gap_extend_percent(2) =  70
  gap_extend_percent(3) =  40
  gap_extend_percent(4) =  85

  go(0:m) = gap_open * real ( gap_open_percent(0:m) ) / 100.0E+00
  ge(0:m) = gap_extend * real ( gap_extend_percent(0:m) ) / 100.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02:'
  write ( *, '(a)' ) '  Forward algorithms using quadratic space:'
  write ( *, '(a)' ) '  PS_GG_FSQ determines an alignment score;'
  write ( *, '(a)' ) '  PS_GG_FPQ determines an alignment path.'

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' I  A  B'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(i3,2x,a1)' ) i, b(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Matching Scores:'
  write ( *, '(a)' ) ' '
  do i = 0, m
    if ( i == 0 ) then
      write ( *, '(i3,13f5.1)' ) i, ( 0.0E+00, j = 0, n )
    else
      write ( *, '(i3,13f5.1)' ) i, 0.0E+00, &
        ( ps_score_simple ( i, b(j) ), j = 1, n )
    end if
  end do

  m1 = 0
  m2 = m
  n1 = 0
  n2 = n
  base = 0.0E+00

  call ps_gg_fsq ( b, m, m1, m2, n, n1, n2, ps_score_simple, go, ge, &
    base, lds, s, e, f, t )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  PS_GG_FSQ:'
  write ( *, '(a)' ) ' '
  write ( *, '(3x,13i5)' ) ( j, j = n1, n2 )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SF:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13f5.1)' ) i, s(i,n1:n2)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'EF:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13f5.1)' ) i, e(i,n1:n2)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FF:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13f5.1)' ) i, f(i,n1:n2)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TF:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13i5)' ) i, t(i,n1:n2)
  end do

  call ps_gg_fpq ( m, m1, m2, n, n1, n2, lds, t, npath, pathi, pathj )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Matching path:'
  write ( *, '(a)' ) ' '
  do i = 1, npath
    write ( *, '(3i6)' ) i, pathi(i), pathj(i)
  end do

  call ps_gg_match_print ( b, m, n, npath, pathi, pathj, ps_score_simple, &
    go, ge )

  return
end
subroutine test03 ( )

!*****************************************************************************80
!
!! TEST03 tests PS_GG_BPQ, PS_GG_BSQ.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 December 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: m = 4
  integer, parameter :: n = 5
  integer, parameter :: lds = m

  character b(n)
  real base
  real e(0:lds,0:n)
  real f(0:lds,0:n)
  real, parameter :: gap_extend = -0.5E+00
  integer gap_extend_percent(0:m)
  real, parameter :: gap_open = -2.0E+00
  integer gap_open_percent(0:m)
  real ge(0:m)
  real go(0:m)
  integer i
  integer j
  integer m1
  integer m2
  integer n1
  integer n2
  integer npath
  integer pathi(m+n+1)
  integer pathj(m+n+1)
  real, external :: ps_score_simple
  real s(0:lds,0:n)
  integer t(0:lds,0:n)

  b(1) = 'A'
  b(2) = 'C'
  b(3) = 'G'
  b(4) = 'C'
  b(5) = 'T'

  gap_open_percent(0) =  95
  gap_open_percent(1) =  90
  gap_open_percent(2) =  75
  gap_open_percent(3) =  50
  gap_open_percent(4) =  25

  gap_extend_percent(0) = 100
  gap_extend_percent(1) =  80
  gap_extend_percent(2) =  70
  gap_extend_percent(3) =  40
  gap_extend_percent(4) =  85

  go(0:m) = gap_open * real ( gap_open_percent(0:m) ) / 100.0E+00
  ge(0:m) = gap_extend * real ( gap_extend_percent(0:m) ) / 100.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03:'
  write ( *, '(a)' ) '  Backward algorithms using quadratic space:'
  write ( *, '(a)' ) '  PS_GG_BSQ determines an alignment score;'
  write ( *, '(a)' ) '  PS_GG_BPQ determines an alignment path.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Verify that these results match those for'
  write ( *, '(a)' ) '  PS_GG_FSQ/PS_GG_FPQ on the same data.'

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' I  A  B'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(i3,2x,a1)' ) i, b(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Matching Scores:'
  write ( *, '(a)' ) ' '
  do i = 0, m
    if ( i == 0 ) then
      write ( *, '(i3,13f5.1)' ) i, ( 0.0E+00, j = 0, n )
    else
      write ( *, '(i3,13f5.1)' ) i, 0.0E+00, &
        ( ps_score_simple ( i, b(j) ), j = 1, n )
    end if
  end do

  m1 = 0
  m2 = m
  n1 = 0
  n2 = n
  base = 0.0E+00

  call ps_gg_bsq ( b, m, m1, m2, n, n1, n2, ps_score_simple, go, ge, &
    base, lds, s, e, f, t )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  PS_GG_BSQ:'
  write ( *, '(a)' ) ' '
  write ( *, '(3x,13i5)' ) ( j, j = n1, n2 )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SB:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13f5.1)' ) i, s(i,n1:n2)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'EB:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13f5.1)' ) i, e(i,n1:n2)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FB:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13f5.1)' ) i, f(i,n1:n2)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TB:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13i5)' ) i, t(i,n1:n2)
  end do

  call ps_gg_bpq ( m, m1, m2, n, n1, n2, lds, t, npath, pathi, pathj )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Matching path:'
  write ( *, '(a)' ) ' '
  do i = 1, npath
    write ( *, '(3i6)' ) i, pathi(i), pathj(i)
  end do

  call ps_gg_match_print ( b, m, n, npath, pathi, pathj, ps_score_simple, &
    go, ge )

  return
end
subroutine test04 ( )

!*****************************************************************************80
!
!! TEST04 tests PS_GG_RPL, PS_GG_FSL.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 December 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: m = 4
  integer, parameter :: n = 5
  integer, parameter :: lds = m

  character b(n)
  real base
  real e(0:lds,0:n)
  real ev(0:n)
  real f(0:lds,0:n)
  real fv(0:n)
  real, parameter :: gap_extend = -0.5E+00
  integer gap_extend_percent(0:m)
  real, parameter :: gap_open = -2.0E+00
  integer gap_open_percent(0:m)
  real ge(0:m)
  real go(0:m)
  integer i
  integer j
  integer m1
  integer m2
  integer n1
  integer n2
  integer npath
  integer pathi(m+n+1)
  integer pathj(m+n+1)
  real, external :: ps_score_simple
  real s(0:lds,0:n)
  real sv(0:n)
  integer t(0:lds,0:n)
  integer tv(0:n)

  b(1) = 'A'
  b(2) = 'C'
  b(3) = 'G'
  b(4) = 'C'
  b(5) = 'T'

  gap_open_percent(0) =  95
  gap_open_percent(1) =  90
  gap_open_percent(2) =  75
  gap_open_percent(3) =  50
  gap_open_percent(4) =  25

  gap_extend_percent(0) = 100
  gap_extend_percent(1) =  80
  gap_extend_percent(2) =  70
  gap_extend_percent(3) =  40
  gap_extend_percent(4) =  85

  go(0:m) = gap_open * real ( gap_open_percent(0:m) ) / 100.0E+00
  ge(0:m) = gap_extend * real ( gap_extend_percent(0:m) ) / 100.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04:'
  write ( *, '(a)' ) '  Forward algorithm using linear space:'
  write ( *, '(a)' ) '  PS_GG_FSL determines an alignment score;'
  write ( *, '(a)' ) '  Recursive algorithm using linear space:'
  write ( *, '(a)' ) '  PS_GG_RPL determines an alignment path.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Verify that these results match those for'
  write ( *, '(a)' ) '  PS_GG_FSQ/PS_GG_FPQ and PS_GG_BSQ/PS_GG_BPQ'
  write ( *, '(a)' ) '  on the same data.'

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' I  A  B'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(i3,2x,a1)' ) i, b(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Matching Scores:'
  write ( *, '(a)' ) ' '
  do i = 0, m
    if ( i == 0 ) then
      write ( *, '(i3,13f5.1)' ) i, ( 0.0E+00, j = 0, n )
    else
      write ( *, '(i3,13f5.1)' ) i, 0.0E+00, &
        ( ps_score_simple ( i, b(j) ), j = 1, n )
    end if
  end do

  m1 = 0
  n1 = 0
  n2 = n
  base = 0.0E+00

  do m2 = m1, m

    call ps_gg_fsl ( b, m, m1, m2, n, n1, n2, ps_score_simple, go, ge, &
      base, sv, ev, fv, tv )

    s(m2,0:n) = sv(0:n)
    e(m2,0:n) = ev(0:n)
    f(m2,0:n) = fv(0:n)
    t(m2,0:n) = tv(0:n)

  end do

  m2 = m

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  PS_GG_FSL:'
  write ( *, '(a)' ) ' '
  write ( *, '(3x,13i5)' ) ( j, j = n1, n2 )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SF:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13f5.1)' ) i, s(i,n1:n2)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'EF:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13f5.1)' ) i, e(i,n1:n2)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FF:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13f5.1)' ) i, f(i,n1:n2)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TF:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13i5)' ) i, t(i,n1:n2)
  end do

  call ps_gg_rpl ( b, m, m1, m2, n, n1, n2, ps_score_simple, go, ge, &
    npath, pathi, pathj )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Matching path:'
  write ( *, '(a)' ) ' '
  do i = 1, npath
    write ( *, '(3i6)' ) i, pathi(i), pathj(i)
  end do

  call ps_gg_match_print ( b, m, n, npath, pathi, pathj, ps_score_simple, &
    go, ge )

  return
end
subroutine test05 ( )

!*****************************************************************************80
!
!! TEST05 tests PS_GG_RPL, PS_GG_BSL.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 December 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: m = 4
  integer, parameter :: n = 5
  integer, parameter :: lds = m

  character b(n)
  real base
  real e(0:lds,0:n)
  real ev(0:n)
  real f(0:lds,0:n)
  real fv(0:n)
  real, parameter :: gap_extend = -0.5E+00
  integer gap_extend_percent(0:m)
  real, parameter :: gap_open = -2.0E+00
  integer gap_open_percent(0:m)
  real ge(0:m)
  real go(0:m)
  integer i
  integer j
  integer m1
  integer m2
  integer n1
  integer n2
  integer npath
  integer pathi(m+n+1)
  integer pathj(m+n+1)
  real, external :: ps_score_simple
  real s(0:lds,0:n)
  real sv(0:n)
  integer t(0:lds,0:n)
  integer tv(0:n)

  b(1) = 'A'
  b(2) = 'C'
  b(3) = 'G'
  b(4) = 'C'
  b(5) = 'T'

  gap_open_percent(0) =  95
  gap_open_percent(1) =  90
  gap_open_percent(2) =  75
  gap_open_percent(3) =  50
  gap_open_percent(4) =  25

  gap_extend_percent(0) = 100
  gap_extend_percent(1) =  80
  gap_extend_percent(2) =  70
  gap_extend_percent(3) =  40
  gap_extend_percent(4) =  85

  go(0:m) = gap_open * real ( gap_open_percent(0:m) ) / 100.0E+00
  ge(0:m) = gap_extend * real ( gap_extend_percent(0:m) ) / 100.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05:'
  write ( *, '(a)' ) '  Backward algorithm using linear space:'
  write ( *, '(a)' ) '  PS_GG_BSL determines an alignment score;'
  write ( *, '(a)' ) '  Recursive algorithm using linear space:'
  write ( *, '(a)' ) '  PS_GG_RPL determines an alignment path.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Verify that these results match those for'
  write ( *, '(a)' ) '  PS_GG_FSQ/PS_GG_FPQ, PS_GG_BSQ/PS_GG_BPQ'
  write ( *, '(a)' ) '  and PS_GG_FSL/PS_GG_RPL on the same data.'

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' I  A  B'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(i3,2x,a1)' ) i, b(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Matching Scores:'
  write ( *, '(a)' ) ' '
  do i = 0, m
    if ( i == 0 ) then
      write ( *, '(i3,13f5.1)' ) i, ( 0.0E+00, j = 0, n )
    else
      write ( *, '(i3,13f5.1)' ) i, 0.0E+00, &
        ( ps_score_simple ( i, b(j) ), j = 1, n )
    end if
  end do

  m1 = 0
  m2 = m
  n1 = 0
  n2 = n
  base = 0.0E+00

  do m1 = 0, m2

    call ps_gg_bsl ( b, m, m1, m2, n, n1, n2, ps_score_simple, go, ge, &
      base, sv, ev, fv, tv )

    s(m1,0:n) = sv(0:n)
    e(m1,0:n) = ev(0:n)
    f(m1,0:n) = fv(0:n)
    t(m1,0:n) = tv(0:n)

  end do

  m1 = 0

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  PS_GG_BSL:'
  write ( *, '(a)' ) ' '
  write ( *, '(3x,13i5)' ) ( j, j = n1, n2 )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SB:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13f5.1)' ) i, s(i,n1:n2)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'EB:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13f5.1)' ) i, e(i,n1:n2)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FB:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13f5.1)' ) i, f(i,n1:n2)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TB:'
  write ( *, '(a)' ) ' '
  do i = m1, m2
    write ( *, '(i3,13i5)' ) i, t(i,n1:n2)
  end do

  call ps_gg_rpl ( b, m, m1, m2, n, n1, n2, ps_score_simple, go, ge, &
    npath, pathi, pathj )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Matching path:'
  write ( *, '(a)' ) ' '
  do i = 1, npath
    write ( *, '(3i6)' ) i, pathi(i), pathj(i)
  end do

  call ps_gg_match_print ( b, m, n, npath, pathi, pathj, ps_score_simple, &
    go, ge )

  return
end
function ps_score_simple ( p1, c2 )

!*****************************************************************************80
!
!! PS_SCORE_SIMPLE computes a single entry profile/sequence matching score.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    31 March 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer P1, the position in the profile.
!
!    Input, character C2, the character in the sequence that is
!    to be matched with profile position P1.
!
!    Output, real PS_SCORE_SIMPLE, the score for matching the given profile
!    position with the sequence character.
!
  implicit none

  character c2
  integer p1
  real ps_score_simple
  real score

  score = 0.0E+00
!
!  Setting P1 = 0 to handle gaps?
!
       if ( p1 == 0 ) then

         if ( c2 == 'A' ) then
      score = 5.0E+00
    else if ( c2 == 'C' ) then
      score = 2.0E+00
    else if ( c2 == 'G' ) then
      score = 4.0E+00
    else if ( c2 == 'T' ) then
      score = 1.0E+00
    else
      score = 0.0E+00
    end if

  else if ( p1 == 1 ) then

         if ( c2 == 'A' ) then
      score = 5.0E+00
    else if ( c2 == 'C' ) then
      score = 2.0E+00
    else if ( c2 == 'G' ) then
      score = 4.0E+00
    else if ( c2 == 'T' ) then
      score = 1.0E+00
    else
      score = 0.0E+00
    end if

  else if ( p1 == 2 ) then

         if ( c2 == 'A' ) then
      score = 2.0E+00
    else if ( c2 == 'C' ) then
      score = 4.0E+00
    else if ( c2 == 'G' ) then
      score = 2.0E+00
    else if ( c2 == 'T' ) then
      score = 1.0E+00
    else
      score = 0.0E+00
    end if

  else if ( p1 == 3 ) then

         if ( c2 == 'A' ) then
      score = 4.0E+00
    else if ( c2 == 'C' ) then
      score = 2.0E+00
    else if ( c2 == 'G' ) then
      score = 4.0E+00
    else if ( c2 == 'T' ) then
      score = 1.0E+00
    else
      score = 0.0E+00
    end if

  else if ( p1 == 4 ) then

         if ( c2 == 'A' ) then
      score = 1.0E+00
    else if ( c2 == 'C' ) then
      score = 1.0E+00
    else if ( c2 == 'G' ) then
      score = 1.0E+00
    else if ( c2 == 'T' ) then
      score = 3.0E+00
    else
      score = 0.0E+00
    end if

  else if ( p1 == 5 ) then

         if ( c2 == 'A' ) then
      score = 0.0E+00
    else if ( c2 == 'C' ) then
      score = 0.0E+00
    else if ( c2 == 'G' ) then
      score = 0.0E+00
    else if ( c2 == 'T' ) then
      score = 0.0E+00
    else
      score = 0.0E+00
    end if

  end if

  ps_score_simple = score

  return
end
subroutine test06 ( )

!*****************************************************************************80
!
!! TEST06 tests PS_GG_FPQ, PS_GG_FSQ.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 December 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: n_max = 200
  integer, parameter :: m_max = 200
  integer, parameter :: lds = m_max

  character b(n_max)
  real base
  real e(0:lds,0:n_max)
  real f(0:lds,0:n_max)
  character ( len = 80 ) file_name
  real, parameter :: gap_extend = -0.5E+00
  integer gap_extend_percent(0:m_max)
  real, parameter :: gap_open = -2.0E+00
  integer gap_open_percent(0:m_max)
  real ge(0:m_max)
  real go(0:m_max)
  integer i
  integer ios
  integer :: iunit = 1
  integer j
  integer m
  integer m1
  integer m2
  integer n
  integer n1
  integer n2
  integer npath
  integer pathi(m_max+n_max+1)
  integer pathj(m_max+n_max+1)
  real, external :: profile_score
  real s(0:lds,0:n_max)
  character ( len = 100 ) sequence
  character ( len = 20 ) sequence_name
  integer t(0:lds,0:n_max)
  integer test

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06:'
  write ( *, '(a)' ) '  Forward algorithms using quadratic space:'
  write ( *, '(a)' ) '  PS_GG_FSQ determines an alignment score;'
  write ( *, '(a)' ) '  PS_GG_FPQ determines an alignment path.'

  file_name = 'profile.txt'

  open ( unit = iunit, file = file_name, status = 'old', iostat = ios )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST06 - Fatal error!'
    write ( *, '(a)' ) '  Could not open the data file.'
    return
  end if

  call profile_score_read2 ( gap_open_percent, gap_extend_percent, iunit, &
    m_max, m )

  close ( unit = iunit )

  go(0:m) = gap_open * real ( gap_open_percent(0:m) ) / 100.0E+00
  ge(0:m) = gap_extend * real ( gap_extend_percent(0:m) ) / 100.0E+00

  do test = 1, 4
!
!  Set the sequence name and values.
!
    if ( test == 1 ) then

      sequence_name = 's00657'

      sequence = 'QECYHGNGQSYRGTYSTTVTGRTCQAWSSMTPHSHSRTPE' // &
                 'YYPNAGLIMNYCRNPDAVAAPYCYTRDPGVRWEYCNLTQC' // &
                 'SD'

    else if ( test == 2 ) then

      sequence_name = 'c61545'

      sequence = 'QDCYHGNGQSYRGTSSTTVTGRKCQSWSSMIPHRHQKTPE' // &
                 'SYPNAGLTMNYCRNPDADKSPWCYTTDPRVRWEFCNLKKC' // &
                 'SE'

    else if ( test == 3 ) then

      sequence_name = 'plhu-k'

      sequence = 'QDCYHGDGQSYRGTSSTTTTGKKCQSWSSMTPHRHQKTPE' // &
                 'NYPNAGLTMNYCRNPDADKGPWCFTTDPSVRWEYCNLKKC' // &
                 'S'

    else if ( test == 4 ) then

      sequence_name = 'Made-Up'

      sequence = 'QDCYHGDGQSCATYRGTSSTTTTGKKCQSWSSMTPHRHQC' // &
                 'ATKTPENYPNAGLTMNYCRNPDADKGPWCFTTDPSVRWEY' // &
                 'CNLKKCS'

    end if
!
!  Convert the sequence string to a vector of characters.
!
    n = len_trim ( sequence )

    call s_to_chvec ( sequence, n, b )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'Sequence name: ' // trim ( sequence_name )
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'Sequence:'
    write ( *, '(a)' ) ' '
    do i = 1, n
      write ( *, '(i5,5x,a1)' ) i, b(i)
    end do

    m1 = 0
    m2 = m
    n1 = 0
    n2 = n
    base = 0.0E+00

    call ps_gg_fsq ( b, m, m1, m2, n, n1, n2, profile_score, go, ge, &
      base, lds, s, e, f, t )

    call ps_gg_fpq ( m, m1, m2, n, n1, n2, lds, t, npath, pathi, pathj )

    call ps_gg_match_print ( b, m, n, npath, pathi, pathj, profile_score, &
      go, ge )

  end do

  return
end
subroutine test07 ( )

!*****************************************************************************80
!
!! TEST07 tests PS_GG_BPQ, PS_GG_BSQ.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 December 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: n_max = 200
  integer, parameter :: m_max = 200
  integer, parameter :: lds = m_max

  character b(n_max)
  real base
  real e(0:lds,0:n_max)
  real f(0:lds,0:n_max)
  character ( len = 80 ) file_name
  real, parameter :: gap_extend = -0.5E+00
  integer gap_extend_percent(0:m_max)
  real, parameter :: gap_open = -2.0E+00
  integer gap_open_percent(0:m_max)
  real ge(0:m_max)
  real go(0:m_max)
  integer i
  integer ios
  integer :: iunit = 1
  integer j
  integer m
  integer m1
  integer m2
  integer n
  integer n1
  integer n2
  integer npath
  integer pathi(m_max+n_max+1)
  integer pathj(m_max+n_max+1)
  real, external :: profile_score
  real s(0:lds,0:n_max)
  character ( len = 100 ) sequence
  character ( len = 20 ) sequence_name
  integer t(0:lds,0:n_max)
  integer test

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07:'
  write ( *, '(a)' ) '  Backward algorithms using quadratic space:'
  write ( *, '(a)' ) '  PS_GG_BSQ determines an alignment score;'
  write ( *, '(a)' ) '  PS_GG_BPQ determines an alignment path.'

  file_name = 'profile.txt'

  open ( unit = iunit, file = file_name, status = 'old', iostat = ios )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST07 - Fatal error!'
    write ( *, '(a)' ) '  Could not open the data file.'
    return
  end if

  call profile_score_read2 ( gap_open_percent, gap_extend_percent, iunit, &
    m_max, m )

  close ( unit = iunit )

  go(0:m) = gap_open * real ( gap_open_percent(0:m) ) / 100.0E+00
  ge(0:m) = gap_extend * real ( gap_extend_percent(0:m) ) / 100.0E+00

  do test = 1, 4
!
!  Set the sequence name and values.
!
    if ( test == 1 ) then

      sequence_name = 's00657'

      sequence = 'QECYHGNGQSYRGTYSTTVTGRTCQAWSSMTPHSHSRTPE' // &
                 'YYPNAGLIMNYCRNPDAVAAPYCYTRDPGVRWEYCNLTQC' // &
                 'SD'

    else if ( test == 2 ) then

      sequence_name = 'c61545'

      sequence = 'QDCYHGNGQSYRGTSSTTVTGRKCQSWSSMIPHRHQKTPE' // &
                 'SYPNAGLTMNYCRNPDADKSPWCYTTDPRVRWEFCNLKKC' // &
                 'SE'

    else if ( test == 3 ) then

      sequence_name = 'plhu-k'

      sequence = 'QDCYHGDGQSYRGTSSTTTTGKKCQSWSSMTPHRHQKTPE' // &
                 'NYPNAGLTMNYCRNPDADKGPWCFTTDPSVRWEYCNLKKC' // &
                 'S'

    else if ( test == 4 ) then

      sequence_name = 'Made-Up'

      sequence = 'QDCYHGDGQSCATYRGTSSTTTTGKKCQSWSSMTPHRHQC' // &
                 'ATKTPENYPNAGLTMNYCRNPDADKGPWCFTTDPSVRWEY' // &
                 'CNLKKCS'

    end if
!
!  Convert the sequence string to a vector of characters.
!
    n = len_trim ( sequence )

    call s_to_chvec ( sequence, n, b )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'Sequence name: ' // trim ( sequence_name )
    write ( *, '(a)' ) ' '

    m1 = 0
    m2 = m
    n1 = 0
    n2 = n
    base = 0.0E+00

    call ps_gg_bsq ( b, m, m1, m2, n, n1, n2, profile_score, go, ge, &
      base, lds, s, e, f, t )

    call ps_gg_bpq ( m, m1, m2, n, n1, n2, lds, t, npath, pathi, pathj )

    call ps_gg_match_print ( b, m, n, npath, pathi, pathj, profile_score, &
      go, ge )

  end do

  return
end
subroutine test08 ( )

!*****************************************************************************80
!
!! TEST08 tests PS_GG_FSL, PS_GG_RPL.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 December 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: n_max = 200
  integer, parameter :: m_max = 200

  character b(n_max)
  character ( len = 80 ) file_name
  real, parameter :: gap_extend = -0.5E+00
  integer gap_extend_percent(0:m_max)
  real, parameter :: gap_open = -2.0E+00
  integer gap_open_percent(0:m_max)
  real ge(0:m_max)
  real go(0:m_max)
  integer i
  integer ios
  integer :: iunit = 1
  integer m
  integer m1
  integer m2
  integer n
  integer n1
  integer n2
  integer npath
  integer pathi(m_max+n_max+1)
  integer pathj(m_max+n_max+1)
  real, external :: profile_score
  character ( len = 100 ) sequence
  character ( len = 20 ) sequence_name
  integer test

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08:'
  write ( *, '(a)' ) '  Forward algorithm using linear space:'
  write ( *, '(a)' ) '  PS_GG_FSL determines an alignment score;'
  write ( *, '(a)' ) '  Recursive algorithm using linear space:'
  write ( *, '(a)' ) '  PS_GG_RPL determines an alignment path.'
  write ( *, '(a)' ) ' '

  file_name = 'profile.txt'

  open ( unit = iunit, file = file_name, status = 'old', iostat = ios )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST08 - Fatal error!'
    write ( *, '(a)' ) '  Could not open the data file.'
    return
  end if

  call profile_score_read2 ( gap_open_percent, gap_extend_percent, iunit, &
    m_max, m )

  close ( unit = iunit )

  go(0:m) = gap_open * real ( gap_open_percent(0:m) ) / 100.0E+00
  ge(0:m) = gap_extend * real ( gap_extend_percent(0:m) ) / 100.0E+00

  do test = 1, 4
!
!  Set the sequence name and values.
!
    if ( test == 1 ) then

      sequence_name = 's00657'

      sequence = 'QECYHGNGQSYRGTYSTTVTGRTCQAWSSMTPHSHSRTPE' // &
                 'YYPNAGLIMNYCRNPDAVAAPYCYTRDPGVRWEYCNLTQC' // &
                 'SD'

    else if ( test == 2 ) then

      sequence_name = 'c61545'

      sequence = 'QDCYHGNGQSYRGTSSTTVTGRKCQSWSSMIPHRHQKTPE' // &
                 'SYPNAGLTMNYCRNPDADKSPWCYTTDPRVRWEFCNLKKC' // &
                 'SE'

    else if ( test == 3 ) then

      sequence_name = 'plhu-k'

      sequence = 'QDCYHGDGQSYRGTSSTTTTGKKCQSWSSMTPHRHQKTPE' // &
                 'NYPNAGLTMNYCRNPDADKGPWCFTTDPSVRWEYCNLKKC' // &
                 'S'

    else if ( test == 4 ) then

      sequence_name = 'Made-Up'

      sequence = 'QDCYHGDGQSCATYRGTSSTTTTGKKCQSWSSMTPHRHQC' // &
                 'ATKTPENYPNAGLTMNYCRNPDADKGPWCFTTDPSVRWEY' // &
                 'CNLKKCS'

    end if
!
!  Convert the sequence string to a vector of characters.
!
    n = len_trim ( sequence )

    call s_to_chvec ( sequence, n, b )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'Sequence name: ' // trim ( sequence_name )
    write ( *, '(a)' ) ' '

    m1 = 0
    m2 = m
    n1 = 0
    n2 = n

    call ps_gg_rpl ( b, m, m1, m2, n, n1, n2, profile_score, go, ge, &
      npath, pathi, pathj )

    call ps_gg_match_print ( b, m, n, npath, pathi, pathj, profile_score, &
      go, ge )

  end do

  return
end
function profile_score ( p1, c2 )

!*****************************************************************************80
!
!! PROFILE_SCORE computes a single entry profile/sequence matching score.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    10 July 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer P1, the position in the profile.
!
!    Input, character C2, the character in the sequence that is
!    to be matched with profile position P1.
!
!    Output, real PROFILE_SCORE, the score for matching the given profile
!    position with the sequence character.
!
  implicit none

  integer, parameter :: acid_num = 23
  logical, parameter :: debug = .false.
  integer, parameter :: iunit = 1
  integer, parameter :: m_max = 200

  integer a_to_i4
  character, save, dimension ( acid_num ) :: acid_code
  integer, save, dimension ( 26 ) :: acid_index
  character c2
  character, save, dimension ( m_max ) :: conserved
  integer, save, dimension ( m_max ) :: entropy
  character ( len = 80 ) file_name
  integer, save, dimension ( m_max ) :: gap_extend_percent
  integer, save, dimension ( m_max ) :: gap_open_percent
  integer i
  character i_to_a
  integer ic2
  integer index_c2
  logical, save :: initialized = .false.
  integer ios
  integer, save :: m
  integer p1
  real profile_score
  integer, save, dimension ( m_max, acid_num ) :: score1
  integer, save, dimension ( acid_num ) :: score2

  if ( .not. initialized ) then

    file_name = 'profile.txt'
!
!  Read the data from the file.
!
    open ( unit = iunit, file = file_name, status = 'old', iostat = ios )

    if ( ios /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'PROFILE_SCORE - Fatal error!'
      write ( *, '(a)' ) '  Could not open the data file.'
      return
    end if

    call profile_score_read ( acid_code, acid_num, conserved, entropy, &
      gap_open_percent, gap_extend_percent, iunit, m_max, m, score1, score2 )

    close ( unit = iunit )

    call a_index ( acid_num, acid_code, acid_index )

    if ( debug ) then

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '    I    Acid_code'
      write ( *, '(a)' ) ' '
      do i = 1, acid_num
        write ( *, '(i5,5x,a1)' ) i, acid_code(i)
      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '    I    Acid_index'
      write ( *, '(a)' ) ' '
      do i = 1, 26
        write ( *, '(i5,5x,a1,4x,i5)' ) i, i_to_a(i), acid_index(i)
      end do

    end if

    initialized = .true.

  end if

  ic2 = a_to_i4 ( c2 )

  index_c2 = acid_index(ic2)

  if ( index_c2 <= 0 ) then
    profile_score = 0.0E+00
  else if ( p1 <= 0 ) then
    profile_score = 0.0E+00
  else
    profile_score = score1(p1,index_c2)
  end if

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  integer values(8)
  integer y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end