subroutine c1f2kb ( ido, l1, na, cc, in1, ch, in2, wa )

!*****************************************************************************80
!
!! C1F2KB is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(in1,l1,ido,2)
  real ch(in2,l1,2,ido)
  real chold1
  real chold2
  integer i
  integer k
  integer na
  real ti2
  real tr2
  real wa(ido,1,2)

  if ( 1 < ido .or. na == 1 ) then

    do k = 1, l1
      ch(1,k,1,1) = cc(1,k,1,1) + cc(1,k,1,2)
      ch(1,k,2,1) = cc(1,k,1,1) - cc(1,k,1,2)
      ch(2,k,1,1) = cc(2,k,1,1) + cc(2,k,1,2)
      ch(2,k,2,1) = cc(2,k,1,1) - cc(2,k,1,2)
    end do

    do i = 2, ido
      do k = 1, l1

        ch(1,k,1,i) = cc(1,k,i,1) + cc(1,k,i,2)
        tr2         = cc(1,k,i,1) - cc(1,k,i,2)
        ch(2,k,1,i) = cc(2,k,i,1) + cc(2,k,i,2)
        ti2         = cc(2,k,i,1) - cc(2,k,i,2)

        ch(2,k,2,i) = wa(i,1,1) * ti2 + wa(i,1,2) * tr2
        ch(1,k,2,i) = wa(i,1,1) * tr2 - wa(i,1,2) * ti2

      end do
    end do

  else

    do k = 1, l1

      chold1      = cc(1,k,1,1) + cc(1,k,1,2)
      cc(1,k,1,2) = cc(1,k,1,1) - cc(1,k,1,2)
      cc(1,k,1,1) = chold1

      chold2      = cc(2,k,1,1) + cc(2,k,1,2)
      cc(2,k,1,2) = cc(2,k,1,1) - cc(2,k,1,2)
      cc(2,k,1,1) = chold2

    end do

  end if

  return
end
subroutine c1f2kf ( ido, l1, na, cc, in1, ch, in2, wa )

!*****************************************************************************80
!
!! C1F2KF is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(in1,l1,ido,2)
  real ch(in2,l1,2,ido)
  real chold1
  real chold2
  integer i
  integer k
  integer na
  real sn
  real ti2
  real tr2
  real wa(ido,1,2)

  if ( 1 < ido ) then

    do k = 1, l1
      ch(1,k,1,1) = cc(1,k,1,1) + cc(1,k,1,2)
      ch(1,k,2,1) = cc(1,k,1,1) - cc(1,k,1,2)
      ch(2,k,1,1) = cc(2,k,1,1) + cc(2,k,1,2)
      ch(2,k,2,1) = cc(2,k,1,1) - cc(2,k,1,2)
    end do

    do i = 2, ido
      do k = 1, l1
        ch(1,k,1,i) = cc(1,k,i,1) + cc(1,k,i,2)
        tr2         = cc(1,k,i,1) - cc(1,k,i,2)
        ch(2,k,1,i) = cc(2,k,i,1) + cc(2,k,i,2)
        ti2         = cc(2,k,i,1) - cc(2,k,i,2)
        ch(2,k,2,i) = wa(i,1,1) * ti2 - wa(i,1,2) * tr2
        ch(1,k,2,i) = wa(i,1,1) * tr2 + wa(i,1,2) * ti2
      end do
    end do

  else if ( na == 1 ) then

    sn = 1.0E+00 / real ( 2 * l1 )

    do k = 1, l1
      ch(1,k,1,1) = sn * ( cc(1,k,1,1) + cc(1,k,1,2) )
      ch(1,k,2,1) = sn * ( cc(1,k,1,1) - cc(1,k,1,2) )
      ch(2,k,1,1) = sn * ( cc(2,k,1,1) + cc(2,k,1,2) )
      ch(2,k,2,1) = sn * ( cc(2,k,1,1) - cc(2,k,1,2) )
    end do

  else

    sn = 1.0E+00 / real ( 2 * l1 )

    do k = 1, l1

      chold1      = sn * ( cc(1,k,1,1) + cc(1,k,1,2) )
      cc(1,k,1,2) = sn * ( cc(1,k,1,1) - cc(1,k,1,2) )
      cc(1,k,1,1) = chold1

      chold2      = sn * ( cc(2,k,1,1) + cc(2,k,1,2) )
      cc(2,k,1,2) = sn * ( cc(2,k,1,1) - cc(2,k,1,2) )
      cc(2,k,1,1) = chold2

    end do

  end if

  return
end
subroutine c1f3kb ( ido, l1, na, cc, in1, ch, in2, wa )

!*****************************************************************************80
!
!! C1F3KB is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(in1,l1,ido,3)
  real ch(in2,l1,3,ido)
  real ci2
  real ci3
  real cr2
  real cr3
  real di2
  real di3
  real dr2
  real dr3
  integer i
  integer k
  integer na
  real, parameter :: taui =  0.866025403784439E+00
  real, parameter :: taur = -0.5E+00
  real ti2
  real tr2
  real wa(ido,2,2)

  if ( 1 < ido .or. na == 1 ) then

    do k = 1, l1

      tr2 = cc(1,k,1,2)+cc(1,k,1,3)
      cr2 = cc(1,k,1,1)+taur*tr2
      ch(1,k,1,1) = cc(1,k,1,1)+tr2
      ti2 = cc(2,k,1,2)+cc(2,k,1,3)
      ci2 = cc(2,k,1,1)+taur*ti2
      ch(2,k,1,1) = cc(2,k,1,1)+ti2
      cr3 = taui*(cc(1,k,1,2)-cc(1,k,1,3))
      ci3 = taui*(cc(2,k,1,2)-cc(2,k,1,3))

      ch(1,k,2,1) = cr2 - ci3
      ch(1,k,3,1) = cr2 + ci3
      ch(2,k,2,1) = ci2 + cr3
      ch(2,k,3,1) = ci2 - cr3

    end do
  
    do i = 2, ido
      do k = 1, l1
        tr2 = cc(1,k,i,2)+cc(1,k,i,3)
        cr2 = cc(1,k,i,1)+taur*tr2
        ch(1,k,1,i) = cc(1,k,i,1)+tr2
        ti2 = cc(2,k,i,2)+cc(2,k,i,3)
        ci2 = cc(2,k,i,1)+taur*ti2
        ch(2,k,1,i) = cc(2,k,i,1)+ti2
        cr3 = taui*(cc(1,k,i,2)-cc(1,k,i,3))
        ci3 = taui*(cc(2,k,i,2)-cc(2,k,i,3))

        dr2 = cr2 - ci3
        dr3 = cr2 + ci3
        di2 = ci2 + cr3
        di3 = ci2 - cr3

        ch(2,k,2,i) = wa(i,1,1) * di2 + wa(i,1,2) * dr2
        ch(1,k,2,i) = wa(i,1,1) * dr2 - wa(i,1,2) * di2
        ch(2,k,3,i) = wa(i,2,1) * di3 + wa(i,2,2) * dr3
        ch(1,k,3,i) = wa(i,2,1) * dr3 - wa(i,2,2) * di3

      end do
    end do

  else

    do k = 1, l1

      tr2 = cc(1,k,1,2)+cc(1,k,1,3)
      cr2 = cc(1,k,1,1)+taur*tr2
      cc(1,k,1,1) = cc(1,k,1,1)+tr2
      ti2 = cc(2,k,1,2)+cc(2,k,1,3)
      ci2 = cc(2,k,1,1)+taur*ti2
      cc(2,k,1,1) = cc(2,k,1,1)+ti2
      cr3 = taui*(cc(1,k,1,2)-cc(1,k,1,3))
      ci3 = taui*(cc(2,k,1,2)-cc(2,k,1,3))

      cc(1,k,1,2) = cr2 - ci3
      cc(1,k,1,3) = cr2 + ci3
      cc(2,k,1,2) = ci2 + cr3
      cc(2,k,1,3) = ci2 - cr3

    end do

  end if

  return
end
subroutine c1f3kf ( ido, l1, na, cc, in1, ch, in2, wa )

!*****************************************************************************80
!
!! C1F3KF is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(in1,l1,ido,3)
  real ch(in2,l1,3,ido)
  real ci2
  real ci3
  real cr2
  real cr3
  real di2
  real di3
  real dr2
  real dr3
  integer i
  integer k
  integer na
  real sn
  real, parameter :: taui = -0.866025403784439E+00
  real, parameter :: taur = -0.5E+00
  real ti2
  real tr2
  real wa(ido,2,2)

  if ( 1 < ido ) then

    do k = 1, l1
      tr2 = cc(1,k,1,2)+cc(1,k,1,3)
      cr2 = cc(1,k,1,1)+taur*tr2
      ch(1,k,1,1) = cc(1,k,1,1)+tr2
      ti2 = cc(2,k,1,2)+cc(2,k,1,3)
      ci2 = cc(2,k,1,1)+taur*ti2
      ch(2,k,1,1) = cc(2,k,1,1)+ti2
      cr3 = taui*(cc(1,k,1,2)-cc(1,k,1,3))
      ci3 = taui*(cc(2,k,1,2)-cc(2,k,1,3))

      ch(1,k,2,1) = cr2 - ci3
      ch(1,k,3,1) = cr2 + ci3
      ch(2,k,2,1) = ci2 + cr3
      ch(2,k,3,1) = ci2 - cr3

    end do

    do i = 2, ido
      do k = 1, l1

        tr2 = cc(1,k,i,2)+cc(1,k,i,3)
        cr2 = cc(1,k,i,1)+taur*tr2
        ch(1,k,1,i) = cc(1,k,i,1)+tr2
        ti2 = cc(2,k,i,2)+cc(2,k,i,3)
        ci2 = cc(2,k,i,1)+taur*ti2
        ch(2,k,1,i) = cc(2,k,i,1)+ti2
        cr3 = taui*(cc(1,k,i,2)-cc(1,k,i,3))
        ci3 = taui*(cc(2,k,i,2)-cc(2,k,i,3))

        dr2 = cr2 - ci3
        dr3 = cr2 + ci3
        di2 = ci2 + cr3
        di3 = ci2 - cr3

        ch(2,k,2,i) = wa(i,1,1) * di2 - wa(i,1,2) * dr2
        ch(1,k,2,i) = wa(i,1,1) * dr2 + wa(i,1,2) * di2
        ch(2,k,3,i) = wa(i,2,1) * di3 - wa(i,2,2) * dr3
        ch(1,k,3,i) = wa(i,2,1) * dr3 + wa(i,2,2) * di3

       end do
    end do

  else if ( na == 1 ) then

    sn = 1.0E+00 / real ( 3 * l1 )

    do k = 1, l1
      tr2 = cc(1,k,1,2)+cc(1,k,1,3)
      cr2 = cc(1,k,1,1)+taur*tr2
      ch(1,k,1,1) = sn*(cc(1,k,1,1)+tr2)
      ti2 = cc(2,k,1,2)+cc(2,k,1,3)
      ci2 = cc(2,k,1,1)+taur*ti2
      ch(2,k,1,1) = sn*(cc(2,k,1,1)+ti2)
      cr3 = taui*(cc(1,k,1,2)-cc(1,k,1,3))
      ci3 = taui*(cc(2,k,1,2)-cc(2,k,1,3))

      ch(1,k,2,1) = sn*(cr2-ci3)
      ch(1,k,3,1) = sn*(cr2+ci3)
      ch(2,k,2,1) = sn*(ci2+cr3)
      ch(2,k,3,1) = sn*(ci2-cr3)

    end do

  else

    sn = 1.0E+00 / real ( 3 * l1 )

    do k = 1, l1
      tr2 = cc(1,k,1,2)+cc(1,k,1,3)
      cr2 = cc(1,k,1,1)+taur*tr2
      cc(1,k,1,1) = sn*(cc(1,k,1,1)+tr2)
      ti2 = cc(2,k,1,2)+cc(2,k,1,3)
      ci2 = cc(2,k,1,1)+taur*ti2
      cc(2,k,1,1) = sn*(cc(2,k,1,1)+ti2)
      cr3 = taui*(cc(1,k,1,2)-cc(1,k,1,3))
      ci3 = taui*(cc(2,k,1,2)-cc(2,k,1,3))

      cc(1,k,1,2) = sn*(cr2-ci3)
      cc(1,k,1,3) = sn*(cr2+ci3)
      cc(2,k,1,2) = sn*(ci2+cr3)
      cc(2,k,1,3) = sn*(ci2-cr3)

    end do

  end if

  return
end
subroutine c1f4kb ( ido, l1, na, cc, in1, ch, in2, wa )

!*****************************************************************************80
!
!! C1F4KB is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(in1,l1,ido,4)
  real ch(in2,l1,4,ido)
  real ci2
  real ci3
  real ci4
  real cr2
  real cr3
  real cr4
  integer i
  integer k
  integer na
  real ti1
  real ti2
  real ti3
  real ti4
  real tr1
  real tr2
  real tr3
  real tr4
  real wa(ido,3,2)

  if ( 1 < ido .or. na == 1 ) then

    do k = 1, l1
      ti1 = cc(2,k,1,1)-cc(2,k,1,3)
      ti2 = cc(2,k,1,1)+cc(2,k,1,3)
      tr4 = cc(2,k,1,4)-cc(2,k,1,2)
      ti3 = cc(2,k,1,2)+cc(2,k,1,4)
      tr1 = cc(1,k,1,1)-cc(1,k,1,3)
      tr2 = cc(1,k,1,1)+cc(1,k,1,3)
      ti4 = cc(1,k,1,2)-cc(1,k,1,4)
      tr3 = cc(1,k,1,2)+cc(1,k,1,4)
      ch(1,k,1,1) = tr2+tr3
      ch(1,k,3,1) = tr2-tr3
      ch(2,k,1,1) = ti2+ti3
      ch(2,k,3,1) = ti2-ti3
      ch(1,k,2,1) = tr1+tr4
      ch(1,k,4,1) = tr1-tr4
      ch(2,k,2,1) = ti1+ti4
      ch(2,k,4,1) = ti1-ti4
    end do

    do i = 2, ido
      do k = 1, l1
        ti1 = cc(2,k,i,1)-cc(2,k,i,3)
        ti2 = cc(2,k,i,1)+cc(2,k,i,3)
        ti3 = cc(2,k,i,2)+cc(2,k,i,4)
        tr4 = cc(2,k,i,4)-cc(2,k,i,2)
        tr1 = cc(1,k,i,1)-cc(1,k,i,3)
        tr2 = cc(1,k,i,1)+cc(1,k,i,3)
        ti4 = cc(1,k,i,2)-cc(1,k,i,4)
        tr3 = cc(1,k,i,2)+cc(1,k,i,4)
        ch(1,k,1,i) = tr2+tr3
        cr3 = tr2-tr3
        ch(2,k,1,i) = ti2+ti3
        ci3 = ti2-ti3
        cr2 = tr1+tr4
        cr4 = tr1-tr4
        ci2 = ti1+ti4
        ci4 = ti1-ti4

        ch(1,k,2,i) = wa(i,1,1)*cr2-wa(i,1,2)*ci2
        ch(2,k,2,i) = wa(i,1,1)*ci2+wa(i,1,2)*cr2
        ch(1,k,3,i) = wa(i,2,1)*cr3-wa(i,2,2)*ci3
        ch(2,k,3,i) = wa(i,2,1)*ci3+wa(i,2,2)*cr3
        ch(1,k,4,i) = wa(i,3,1)*cr4-wa(i,3,2)*ci4
        ch(2,k,4,i) = wa(i,3,1)*ci4+wa(i,3,2)*cr4

       end do
    end do

  else

    do k = 1, l1
       ti1 = cc(2,k,1,1)-cc(2,k,1,3)
       ti2 = cc(2,k,1,1)+cc(2,k,1,3)
       tr4 = cc(2,k,1,4)-cc(2,k,1,2)
       ti3 = cc(2,k,1,2)+cc(2,k,1,4)
       tr1 = cc(1,k,1,1)-cc(1,k,1,3)
       tr2 = cc(1,k,1,1)+cc(1,k,1,3)
       ti4 = cc(1,k,1,2)-cc(1,k,1,4)
       tr3 = cc(1,k,1,2)+cc(1,k,1,4)
       cc(1,k,1,1) = tr2+tr3
       cc(1,k,1,3) = tr2-tr3
       cc(2,k,1,1) = ti2+ti3
       cc(2,k,1,3) = ti2-ti3
       cc(1,k,1,2) = tr1+tr4
       cc(1,k,1,4) = tr1-tr4
       cc(2,k,1,2) = ti1+ti4
       cc(2,k,1,4) = ti1-ti4
    end do

  end if

  return
end
subroutine c1f4kf ( ido, l1, na, cc, in1, ch, in2, wa )

!*****************************************************************************80
!
!! C1F4KF is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(in1,l1,ido,4)
  real ch(in2,l1,4,ido)
  real ci2
  real ci3
  real ci4
  real cr2
  real cr3
  real cr4
  integer i
  integer k
  integer na
  real sn
  real ti1
  real ti2
  real ti3
  real ti4
  real tr1
  real tr2
  real tr3
  real tr4
  real wa(ido,3,2)

  if ( 1 < ido ) then

    do k = 1, l1

      ti1 = cc(2,k,1,1)-cc(2,k,1,3)
      ti2 = cc(2,k,1,1)+cc(2,k,1,3)
      tr4 = cc(2,k,1,2)-cc(2,k,1,4)
      ti3 = cc(2,k,1,2)+cc(2,k,1,4)
      tr1 = cc(1,k,1,1)-cc(1,k,1,3)
      tr2 = cc(1,k,1,1)+cc(1,k,1,3)
      ti4 = cc(1,k,1,4)-cc(1,k,1,2)
      tr3 = cc(1,k,1,2)+cc(1,k,1,4)

      ch(1,k,1,1) = tr2 + tr3
      ch(1,k,3,1) = tr2 - tr3
      ch(2,k,1,1) = ti2 + ti3
      ch(2,k,3,1) = ti2 - ti3
      ch(1,k,2,1) = tr1 + tr4
      ch(1,k,4,1) = tr1 - tr4
      ch(2,k,2,1) = ti1 + ti4
      ch(2,k,4,1) = ti1 - ti4

    end do

    do i = 2, ido
      do k = 1, l1
        ti1 = cc(2,k,i,1)-cc(2,k,i,3)
        ti2 = cc(2,k,i,1)+cc(2,k,i,3)
        ti3 = cc(2,k,i,2)+cc(2,k,i,4)
        tr4 = cc(2,k,i,2)-cc(2,k,i,4)
        tr1 = cc(1,k,i,1)-cc(1,k,i,3)
        tr2 = cc(1,k,i,1)+cc(1,k,i,3)
        ti4 = cc(1,k,i,4)-cc(1,k,i,2)
        tr3 = cc(1,k,i,2)+cc(1,k,i,4)
        ch(1,k,1,i) = tr2+tr3
        cr3 = tr2-tr3
        ch(2,k,1,i) = ti2+ti3
        ci3 = ti2-ti3
        cr2 = tr1+tr4
        cr4 = tr1-tr4
        ci2 = ti1+ti4
        ci4 = ti1-ti4
        ch(1,k,2,i) = wa(i,1,1)*cr2+wa(i,1,2)*ci2
        ch(2,k,2,i) = wa(i,1,1)*ci2-wa(i,1,2)*cr2
        ch(1,k,3,i) = wa(i,2,1)*cr3+wa(i,2,2)*ci3
        ch(2,k,3,i) = wa(i,2,1)*ci3-wa(i,2,2)*cr3
        ch(1,k,4,i) = wa(i,3,1)*cr4+wa(i,3,2)*ci4
        ch(2,k,4,i) = wa(i,3,1)*ci4-wa(i,3,2)*cr4
      end do
    end do

  else if ( na == 1 ) then

    sn = 1.0E+00 / real ( 4 * l1 )

    do k = 1, l1
      ti1 = cc(2,k,1,1)-cc(2,k,1,3)
      ti2 = cc(2,k,1,1)+cc(2,k,1,3)
      tr4 = cc(2,k,1,2)-cc(2,k,1,4)
      ti3 = cc(2,k,1,2)+cc(2,k,1,4)
      tr1 = cc(1,k,1,1)-cc(1,k,1,3)
      tr2 = cc(1,k,1,1)+cc(1,k,1,3)
      ti4 = cc(1,k,1,4)-cc(1,k,1,2)
      tr3 = cc(1,k,1,2)+cc(1,k,1,4)
      ch(1,k,1,1) = sn*(tr2+tr3)
      ch(1,k,3,1) = sn*(tr2-tr3)
      ch(2,k,1,1) = sn*(ti2+ti3)
      ch(2,k,3,1) = sn*(ti2-ti3)
      ch(1,k,2,1) = sn*(tr1+tr4)
      ch(1,k,4,1) = sn*(tr1-tr4)
      ch(2,k,2,1) = sn*(ti1+ti4)
      ch(2,k,4,1) = sn*(ti1-ti4)
    end do

  else

    sn = 1.0E+00 / real ( 4 * l1 )

    do k = 1, l1
      ti1 = cc(2,k,1,1)-cc(2,k,1,3)
      ti2 = cc(2,k,1,1)+cc(2,k,1,3)
      tr4 = cc(2,k,1,2)-cc(2,k,1,4)
      ti3 = cc(2,k,1,2)+cc(2,k,1,4)
      tr1 = cc(1,k,1,1)-cc(1,k,1,3)
      tr2 = cc(1,k,1,1)+cc(1,k,1,3)
      ti4 = cc(1,k,1,4)-cc(1,k,1,2)
      tr3 = cc(1,k,1,2)+cc(1,k,1,4)
      cc(1,k,1,1) = sn*(tr2+tr3)
      cc(1,k,1,3) = sn*(tr2-tr3)
      cc(2,k,1,1) = sn*(ti2+ti3)
      cc(2,k,1,3) = sn*(ti2-ti3)
      cc(1,k,1,2) = sn*(tr1+tr4)
      cc(1,k,1,4) = sn*(tr1-tr4)
      cc(2,k,1,2) = sn*(ti1+ti4)
      cc(2,k,1,4) = sn*(ti1-ti4)
    end do

  end if

  return
end
subroutine c1f5kb ( ido, l1, na, cc, in1, ch, in2, wa )

!*****************************************************************************80
!
!! C1F5KB is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(in1,l1,ido,5)
  real ch(in2,l1,5,ido)
  real chold1
  real chold2
  real ci2
  real ci3
  real ci4
  real ci5
  real cr2
  real cr3
  real cr4
  real cr5
  real di2
  real di3
  real di4
  real di5
  real dr2
  real dr3
  real dr4
  real dr5
  integer i
  integer k
  integer na
  real ti2
  real ti3
  real ti4
  real ti5
  real, parameter :: ti11 =  0.9510565162951536E+00
  real, parameter :: ti12 =  0.5877852522924731E+00
  real tr2
  real tr3
  real tr4
  real tr5
  real, parameter :: tr11 =  0.3090169943749474E+00
  real, parameter :: tr12 = -0.8090169943749474E+00
  real wa(ido,4,2)

  if ( 1 < ido .or. na == 1 ) then

    do k = 1, l1
      ti5 = cc(2,k,1,2)-cc(2,k,1,5)
      ti2 = cc(2,k,1,2)+cc(2,k,1,5)
      ti4 = cc(2,k,1,3)-cc(2,k,1,4)
      ti3 = cc(2,k,1,3)+cc(2,k,1,4)
      tr5 = cc(1,k,1,2)-cc(1,k,1,5)
      tr2 = cc(1,k,1,2)+cc(1,k,1,5)
      tr4 = cc(1,k,1,3)-cc(1,k,1,4)
      tr3 = cc(1,k,1,3)+cc(1,k,1,4)
      ch(1,k,1,1) = cc(1,k,1,1)+tr2+tr3
      ch(2,k,1,1) = cc(2,k,1,1)+ti2+ti3
      cr2 = cc(1,k,1,1)+tr11*tr2+tr12*tr3
      ci2 = cc(2,k,1,1)+tr11*ti2+tr12*ti3
      cr3 = cc(1,k,1,1)+tr12*tr2+tr11*tr3
      ci3 = cc(2,k,1,1)+tr12*ti2+tr11*ti3
      cr5 = ti11*tr5+ti12*tr4
      ci5 = ti11*ti5+ti12*ti4
      cr4 = ti12*tr5-ti11*tr4
      ci4 = ti12*ti5-ti11*ti4
      ch(1,k,2,1) = cr2-ci5
      ch(1,k,5,1) = cr2+ci5
      ch(2,k,2,1) = ci2+cr5
      ch(2,k,3,1) = ci3+cr4
      ch(1,k,3,1) = cr3-ci4
      ch(1,k,4,1) = cr3+ci4
      ch(2,k,4,1) = ci3-cr4
      ch(2,k,5,1) = ci2-cr5
    end do

    do i = 2, ido
      do k = 1, l1
        ti5 = cc(2,k,i,2)-cc(2,k,i,5)
        ti2 = cc(2,k,i,2)+cc(2,k,i,5)
        ti4 = cc(2,k,i,3)-cc(2,k,i,4)
        ti3 = cc(2,k,i,3)+cc(2,k,i,4)
        tr5 = cc(1,k,i,2)-cc(1,k,i,5)
        tr2 = cc(1,k,i,2)+cc(1,k,i,5)
        tr4 = cc(1,k,i,3)-cc(1,k,i,4)
        tr3 = cc(1,k,i,3)+cc(1,k,i,4)
        ch(1,k,1,i) = cc(1,k,i,1)+tr2+tr3
        ch(2,k,1,i) = cc(2,k,i,1)+ti2+ti3
        cr2 = cc(1,k,i,1)+tr11*tr2+tr12*tr3
        ci2 = cc(2,k,i,1)+tr11*ti2+tr12*ti3
        cr3 = cc(1,k,i,1)+tr12*tr2+tr11*tr3
        ci3 = cc(2,k,i,1)+tr12*ti2+tr11*ti3
        cr5 = ti11*tr5+ti12*tr4
        ci5 = ti11*ti5+ti12*ti4
        cr4 = ti12*tr5-ti11*tr4
        ci4 = ti12*ti5-ti11*ti4
        dr3 = cr3-ci4
        dr4 = cr3+ci4
        di3 = ci3+cr4
        di4 = ci3-cr4
        dr5 = cr2+ci5
        dr2 = cr2-ci5
        di5 = ci2-cr5
        di2 = ci2+cr5
        ch(1,k,2,i) = wa(i,1,1)*dr2-wa(i,1,2)*di2
        ch(2,k,2,i) = wa(i,1,1)*di2+wa(i,1,2)*dr2
        ch(1,k,3,i) = wa(i,2,1)*dr3-wa(i,2,2)*di3
        ch(2,k,3,i) = wa(i,2,1)*di3+wa(i,2,2)*dr3
        ch(1,k,4,i) = wa(i,3,1)*dr4-wa(i,3,2)*di4
        ch(2,k,4,i) = wa(i,3,1)*di4+wa(i,3,2)*dr4
        ch(1,k,5,i) = wa(i,4,1)*dr5-wa(i,4,2)*di5
        ch(2,k,5,i) = wa(i,4,1)*di5+wa(i,4,2)*dr5
      end do
    end do

  else

    do k = 1, l1
      ti5 = cc(2,k,1,2)-cc(2,k,1,5)
      ti2 = cc(2,k,1,2)+cc(2,k,1,5)
      ti4 = cc(2,k,1,3)-cc(2,k,1,4)
      ti3 = cc(2,k,1,3)+cc(2,k,1,4)
      tr5 = cc(1,k,1,2)-cc(1,k,1,5)
      tr2 = cc(1,k,1,2)+cc(1,k,1,5)
      tr4 = cc(1,k,1,3)-cc(1,k,1,4)
      tr3 = cc(1,k,1,3)+cc(1,k,1,4)
      chold1 = cc(1,k,1,1)+tr2+tr3
      chold2 = cc(2,k,1,1)+ti2+ti3
      cr2 = cc(1,k,1,1)+tr11*tr2+tr12*tr3
      ci2 = cc(2,k,1,1)+tr11*ti2+tr12*ti3
      cr3 = cc(1,k,1,1)+tr12*tr2+tr11*tr3
      ci3 = cc(2,k,1,1)+tr12*ti2+tr11*ti3
      cc(1,k,1,1) = chold1
      cc(2,k,1,1) = chold2
      cr5 = ti11*tr5+ti12*tr4
      ci5 = ti11*ti5+ti12*ti4
      cr4 = ti12*tr5-ti11*tr4
      ci4 = ti12*ti5-ti11*ti4
      cc(1,k,1,2) = cr2-ci5
      cc(1,k,1,5) = cr2+ci5
      cc(2,k,1,2) = ci2+cr5
      cc(2,k,1,3) = ci3+cr4
      cc(1,k,1,3) = cr3-ci4
      cc(1,k,1,4) = cr3+ci4
      cc(2,k,1,4) = ci3-cr4
      cc(2,k,1,5) = ci2-cr5
    end do

  end if

  return
end
subroutine c1f5kf ( ido, l1, na, cc, in1, ch, in2, wa )

!*****************************************************************************80
!
!! C1F5KF is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(in1,l1,ido,5)
  real ch(in2,l1,5,ido)
  real chold1
  real chold2
  real ci2
  real ci3
  real ci4
  real ci5
  real cr2
  real cr3
  real cr4
  real cr5
  real di2
  real di3
  real di4
  real di5
  real dr2
  real dr3
  real dr4
  real dr5
  integer i
  integer k
  integer na
  real sn
  real ti2
  real ti3
  real ti4
  real ti5
  real, parameter :: ti11 = -0.9510565162951536E+00
  real, parameter :: ti12 = -0.5877852522924731E+00
  real tr2
  real tr3
  real tr4
  real tr5
  real, parameter :: tr11 =  0.3090169943749474E+00
  real, parameter :: tr12 = -0.8090169943749474E+00
  real wa(ido,4,2)

  if ( 1 < ido ) then

    do k = 1, l1

      ti5 = cc(2,k,1,2)-cc(2,k,1,5)
      ti2 = cc(2,k,1,2)+cc(2,k,1,5)
      ti4 = cc(2,k,1,3)-cc(2,k,1,4)
      ti3 = cc(2,k,1,3)+cc(2,k,1,4)
      tr5 = cc(1,k,1,2)-cc(1,k,1,5)
      tr2 = cc(1,k,1,2)+cc(1,k,1,5)
      tr4 = cc(1,k,1,3)-cc(1,k,1,4)
      tr3 = cc(1,k,1,3)+cc(1,k,1,4)

      ch(1,k,1,1) = cc(1,k,1,1)+tr2+tr3
      ch(2,k,1,1) = cc(2,k,1,1)+ti2+ti3
      cr2 = cc(1,k,1,1)+tr11*tr2+tr12*tr3
      ci2 = cc(2,k,1,1)+tr11*ti2+tr12*ti3
      cr3 = cc(1,k,1,1)+tr12*tr2+tr11*tr3
      ci3 = cc(2,k,1,1)+tr12*ti2+tr11*ti3
      cr5 = ti11*tr5+ti12*tr4
      ci5 = ti11*ti5+ti12*ti4
      cr4 = ti12*tr5-ti11*tr4
      ci4 = ti12*ti5-ti11*ti4
      ch(1,k,2,1) = cr2-ci5
      ch(1,k,5,1) = cr2+ci5
      ch(2,k,2,1) = ci2+cr5
      ch(2,k,3,1) = ci3+cr4
      ch(1,k,3,1) = cr3-ci4
      ch(1,k,4,1) = cr3+ci4
      ch(2,k,4,1) = ci3-cr4
      ch(2,k,5,1) = ci2-cr5
    end do

    do i = 2, ido
      do k = 1, l1

        ti5 = cc(2,k,i,2)-cc(2,k,i,5)
        ti2 = cc(2,k,i,2)+cc(2,k,i,5)
        ti4 = cc(2,k,i,3)-cc(2,k,i,4)
        ti3 = cc(2,k,i,3)+cc(2,k,i,4)
        tr5 = cc(1,k,i,2)-cc(1,k,i,5)
        tr2 = cc(1,k,i,2)+cc(1,k,i,5)
        tr4 = cc(1,k,i,3)-cc(1,k,i,4)
        tr3 = cc(1,k,i,3)+cc(1,k,i,4)

        ch(1,k,1,i) = cc(1,k,i,1)+tr2+tr3
        ch(2,k,1,i) = cc(2,k,i,1)+ti2+ti3
        cr2 = cc(1,k,i,1)+tr11*tr2+tr12*tr3
        ci2 = cc(2,k,i,1)+tr11*ti2+tr12*ti3
        cr3 = cc(1,k,i,1)+tr12*tr2+tr11*tr3
        ci3 = cc(2,k,i,1)+tr12*ti2+tr11*ti3
        cr5 = ti11*tr5+ti12*tr4
        ci5 = ti11*ti5+ti12*ti4
        cr4 = ti12*tr5-ti11*tr4
        ci4 = ti12*ti5-ti11*ti4
        dr3 = cr3-ci4
        dr4 = cr3+ci4
        di3 = ci3+cr4
        di4 = ci3-cr4
        dr5 = cr2+ci5
        dr2 = cr2-ci5
        di5 = ci2-cr5
        di2 = ci2+cr5
        ch(1,k,2,i) = wa(i,1,1)*dr2+wa(i,1,2)*di2
        ch(2,k,2,i) = wa(i,1,1)*di2-wa(i,1,2)*dr2
        ch(1,k,3,i) = wa(i,2,1)*dr3+wa(i,2,2)*di3
        ch(2,k,3,i) = wa(i,2,1)*di3-wa(i,2,2)*dr3
        ch(1,k,4,i) = wa(i,3,1)*dr4+wa(i,3,2)*di4
        ch(2,k,4,i) = wa(i,3,1)*di4-wa(i,3,2)*dr4
        ch(1,k,5,i) = wa(i,4,1)*dr5+wa(i,4,2)*di5
        ch(2,k,5,i) = wa(i,4,1)*di5-wa(i,4,2)*dr5
      end do
    end do

  else if ( na == 1 ) then

    sn = 1.0E+00 / real ( 5 * l1 )

    do k = 1, l1

      ti5 = cc(2,k,1,2)-cc(2,k,1,5)
      ti2 = cc(2,k,1,2)+cc(2,k,1,5)
      ti4 = cc(2,k,1,3)-cc(2,k,1,4)
      ti3 = cc(2,k,1,3)+cc(2,k,1,4)
      tr5 = cc(1,k,1,2)-cc(1,k,1,5)
      tr2 = cc(1,k,1,2)+cc(1,k,1,5)
      tr4 = cc(1,k,1,3)-cc(1,k,1,4)
      tr3 = cc(1,k,1,3)+cc(1,k,1,4)

      ch(1,k,1,1) = sn*(cc(1,k,1,1)+tr2+tr3)
      ch(2,k,1,1) = sn*(cc(2,k,1,1)+ti2+ti3)

      cr2 = cc(1,k,1,1)+tr11*tr2+tr12*tr3
      ci2 = cc(2,k,1,1)+tr11*ti2+tr12*ti3
      cr3 = cc(1,k,1,1)+tr12*tr2+tr11*tr3
      ci3 = cc(2,k,1,1)+tr12*ti2+tr11*ti3
      cr5 = ti11*tr5+ti12*tr4
      ci5 = ti11*ti5+ti12*ti4
      cr4 = ti12*tr5-ti11*tr4
      ci4 = ti12*ti5-ti11*ti4

      ch(1,k,2,1) = sn*(cr2-ci5)
      ch(1,k,5,1) = sn*(cr2+ci5)
      ch(2,k,2,1) = sn*(ci2+cr5)
      ch(2,k,3,1) = sn*(ci3+cr4)
      ch(1,k,3,1) = sn*(cr3-ci4)
      ch(1,k,4,1) = sn*(cr3+ci4)
      ch(2,k,4,1) = sn*(ci3-cr4)
      ch(2,k,5,1) = sn*(ci2-cr5)

    end do

  else

    sn = 1.0E+00 / real ( 5 * l1 )

    do k = 1, l1

      ti5 = cc(2,k,1,2)-cc(2,k,1,5)
      ti2 = cc(2,k,1,2)+cc(2,k,1,5)
      ti4 = cc(2,k,1,3)-cc(2,k,1,4)
      ti3 = cc(2,k,1,3)+cc(2,k,1,4)
      tr5 = cc(1,k,1,2)-cc(1,k,1,5)
      tr2 = cc(1,k,1,2)+cc(1,k,1,5)
      tr4 = cc(1,k,1,3)-cc(1,k,1,4)
      tr3 = cc(1,k,1,3)+cc(1,k,1,4)

      chold1 = sn*(cc(1,k,1,1)+tr2+tr3)
      chold2 = sn*(cc(2,k,1,1)+ti2+ti3)

      cr2 = cc(1,k,1,1)+tr11*tr2+tr12*tr3
      ci2 = cc(2,k,1,1)+tr11*ti2+tr12*ti3
      cr3 = cc(1,k,1,1)+tr12*tr2+tr11*tr3
      ci3 = cc(2,k,1,1)+tr12*ti2+tr11*ti3

      cc(1,k,1,1) = chold1
      cc(2,k,1,1) = chold2

      cr5 = ti11*tr5+ti12*tr4
      ci5 = ti11*ti5+ti12*ti4
      cr4 = ti12*tr5-ti11*tr4
      ci4 = ti12*ti5-ti11*ti4

      cc(1,k,1,2) = sn*(cr2-ci5)
      cc(1,k,1,5) = sn*(cr2+ci5)
      cc(2,k,1,2) = sn*(ci2+cr5)
      cc(2,k,1,3) = sn*(ci3+cr4)
      cc(1,k,1,3) = sn*(cr3-ci4)
      cc(1,k,1,4) = sn*(cr3+ci4)
      cc(2,k,1,4) = sn*(ci3-cr4)
      cc(2,k,1,5) = sn*(ci2-cr5)

    end do

  end if

  return
end
subroutine c1fgkb ( ido, ip, l1, lid, na, cc, cc1, in1, ch, ch1, in2, wa )

!*****************************************************************************80
!
!! C1FGKB is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer ip
  integer l1
  integer lid

  real cc(in1,l1,ip,ido)
  real cc1(in1,lid,ip)
  real ch(in2,l1,ido,ip)
  real ch1(in2,lid,ip)
  real chold1
  real chold2
  integer i
  integer idlj
  integer ipp2
  integer ipph
  integer j
  integer jc
  integer k
  integer ki
  integer l
  integer lc
  integer na
  real wa(ido,ip-1,2)
  real wai
  real war

  ipp2 = ip + 2
  ipph = ( ip + 1 ) / 2
  do ki = 1, lid
    ch1(1,ki,1) = cc1(1,ki,1)
    ch1(2,ki,1) = cc1(2,ki,1)
  end do

  do j = 2, ipph
    jc = ipp2 - j
    do ki = 1, lid
      ch1(1,ki,j) =  cc1(1,ki,j) + cc1(1,ki,jc)
      ch1(1,ki,jc) = cc1(1,ki,j) - cc1(1,ki,jc)
      ch1(2,ki,j) =  cc1(2,ki,j) + cc1(2,ki,jc)
      ch1(2,ki,jc) = cc1(2,ki,j) - cc1(2,ki,jc)
    end do
  end do

  do j = 2, ipph
    do ki = 1, lid
      cc1(1,ki,1) = cc1(1,ki,1)+ch1(1,ki,j)
      cc1(2,ki,1) = cc1(2,ki,1)+ch1(2,ki,j)
    end do
  end do

  do l = 2, ipph

     lc = ipp2 - l
     do ki = 1, lid
       cc1(1,ki,l) = ch1(1,ki,1)+wa(1,l-1,1)*ch1(1,ki,2)
       cc1(1,ki,lc) = wa(1,l-1,2)*ch1(1,ki,ip)
       cc1(2,ki,l) = ch1(2,ki,1)+wa(1,l-1,1)*ch1(2,ki,2)
       cc1(2,ki,lc) = wa(1,l-1,2)*ch1(2,ki,ip)
     end do

     do j = 3, ipph
       jc = ipp2 - j
       idlj = mod ( ( l - 1 ) * ( j - 1 ), ip )
       war = wa(1,idlj,1)
       wai = wa(1,idlj,2)
       do ki = 1, lid
         cc1(1,ki,l) = cc1(1,ki,l)+war*ch1(1,ki,j)
         cc1(1,ki,lc) = cc1(1,ki,lc)+wai*ch1(1,ki,jc)
         cc1(2,ki,l) = cc1(2,ki,l)+war*ch1(2,ki,j)
         cc1(2,ki,lc) = cc1(2,ki,lc)+wai*ch1(2,ki,jc)
       end do
     end do

  end do

  if ( 1 < ido .or. na == 1 ) then

    do ki = 1, lid
      ch1(1,ki,1) = cc1(1,ki,1)
      ch1(2,ki,1) = cc1(2,ki,1)
    end do

    do j = 2, ipph
      jc = ipp2 - j
      do ki = 1, lid
        ch1(1,ki,j) = cc1(1,ki,j)-cc1(2,ki,jc)
        ch1(1,ki,jc) = cc1(1,ki,j)+cc1(2,ki,jc)
        ch1(2,ki,jc) = cc1(2,ki,j)-cc1(1,ki,jc)
        ch1(2,ki,j) = cc1(2,ki,j)+cc1(1,ki,jc)
      end do
    end do

    if ( ido == 1 ) then
      return
    end if

    do i = 1, ido
      do k = 1, l1
        cc(1,k,1,i) = ch(1,k,i,1)
        cc(2,k,1,i) = ch(2,k,i,1)
      end do
    end do

    do j = 2, ip
      do k = 1, l1
        cc(1,k,j,1) = ch(1,k,1,j)
        cc(2,k,j,1) = ch(2,k,1,j)
      end do
    end do

    do j = 2, ip
      do i = 2, ido
        do k = 1, l1
          cc(1,k,j,i) = wa(i,j-1,1)*ch(1,k,i,j) &
                       -wa(i,j-1,2)*ch(2,k,i,j)
          cc(2,k,j,i) = wa(i,j-1,1)*ch(2,k,i,j) &
                       +wa(i,j-1,2)*ch(1,k,i,j)
        end do
      end do
    end do

  else

    do j = 2, ipph
      jc = ipp2 - j
      do ki = 1, lid
        chold1 = cc1(1,ki,j)-cc1(2,ki,jc)
        chold2 = cc1(1,ki,j)+cc1(2,ki,jc)
        cc1(1,ki,j) = chold1
        cc1(2,ki,jc) = cc1(2,ki,j)-cc1(1,ki,jc)
        cc1(2,ki,j) = cc1(2,ki,j)+cc1(1,ki,jc)
        cc1(1,ki,jc) = chold2
      end do
    end do

  end if

  return
end
subroutine c1fgkf ( ido, ip, l1, lid, na, cc, cc1, in1, ch, ch1, in2, wa )

!*****************************************************************************80
!
!! C1FGKF is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer ip
  integer l1
  integer lid

  real cc(in1,l1,ip,ido)
  real cc1(in1,lid,ip)
  real ch(in2,l1,ido,ip)
  real ch1(in2,lid,ip)
  real chold1
  real chold2
  integer i
  integer idlj
  integer ipp2
  integer ipph
  integer j
  integer jc
  integer k
  integer ki
  integer l
  integer lc
  integer na
  real sn
  real wa(ido,ip-1,2)
  real wai
  real war

  ipp2 = ip+2
  ipph = (ip+1)/2

  do ki = 1, lid
    ch1(1,ki,1) = cc1(1,ki,1)
    ch1(2,ki,1) = cc1(2,ki,1)
  end do

  do j = 2, ipph
    jc = ipp2 - j
    do ki = 1, lid
      ch1(1,ki,j) =  cc1(1,ki,j)+cc1(1,ki,jc)
      ch1(1,ki,jc) = cc1(1,ki,j)-cc1(1,ki,jc)
      ch1(2,ki,j) =  cc1(2,ki,j)+cc1(2,ki,jc)
      ch1(2,ki,jc) = cc1(2,ki,j)-cc1(2,ki,jc)
    end do
  end do

  do j = 2, ipph
    do ki = 1, lid
      cc1(1,ki,1) = cc1(1,ki,1) + ch1(1,ki,j)
      cc1(2,ki,1) = cc1(2,ki,1) + ch1(2,ki,j)
    end do
  end do

  do l = 2, ipph

    lc = ipp2 - l

    do ki = 1, lid
      cc1(1,ki,l)  = ch1(1,ki,1) + wa(1,l-1,1) * ch1(1,ki,2)
      cc1(1,ki,lc) =             - wa(1,l-1,2) * ch1(1,ki,ip)
      cc1(2,ki,l)  = ch1(2,ki,1) + wa(1,l-1,1) * ch1(2,ki,2)
      cc1(2,ki,lc) =             - wa(1,l-1,2) * ch1(2,ki,ip)
    end do

    do j = 3, ipph

      jc = ipp2 - j
      idlj = mod ( ( l - 1 ) * ( j - 1 ), ip )
      war = wa(1,idlj,1)
      wai = -wa(1,idlj,2)

      do ki = 1, lid
        cc1(1,ki,l) = cc1(1,ki,l)+war*ch1(1,ki,j)
        cc1(1,ki,lc) = cc1(1,ki,lc)+wai*ch1(1,ki,jc)
        cc1(2,ki,l) = cc1(2,ki,l)+war*ch1(2,ki,j)
        cc1(2,ki,lc) = cc1(2,ki,lc)+wai*ch1(2,ki,jc)
      end do

    end do

  end do

  if ( 1 < ido ) then

    do ki = 1, lid
      ch1(1,ki,1) = cc1(1,ki,1)
      ch1(2,ki,1) = cc1(2,ki,1)
    end do

    do j = 2, ipph
      jc = ipp2 - j
      do ki = 1, lid
        ch1(1,ki,j) = cc1(1,ki,j)-cc1(2,ki,jc)
        ch1(2,ki,j) = cc1(2,ki,j)+cc1(1,ki,jc)
        ch1(1,ki,jc) = cc1(1,ki,j)+cc1(2,ki,jc)
        ch1(2,ki,jc) = cc1(2,ki,j)-cc1(1,ki,jc)
      end do
    end do

    do i = 1, ido
      do k = 1, l1
        cc(1,k,1,i) = ch(1,k,i,1)
        cc(2,k,1,i) = ch(2,k,i,1)
      end do
    end do

    do j = 2, ip
      do k = 1, l1
        cc(1,k,j,1) = ch(1,k,1,j)
        cc(2,k,j,1) = ch(2,k,1,j)
      end do
    end do

    do j = 2, ip
      do i = 2, ido
        do k = 1, l1
          cc(1,k,j,i) = wa(i,j-1,1)*ch(1,k,i,j) + wa(i,j-1,2)*ch(2,k,i,j)
          cc(2,k,j,i) = wa(i,j-1,1)*ch(2,k,i,j) - wa(i,j-1,2)*ch(1,k,i,j)
        end do
      end do
    end do

  else if ( na == 1 ) then

    sn = 1.0E+00 / real ( ip * l1 )

    do ki = 1, lid
      ch1(1,ki,1) = sn * cc1(1,ki,1)
      ch1(2,ki,1) = sn * cc1(2,ki,1)
    end do

    do j = 2, ipph
      jc = ipp2 - j
      do ki = 1, lid
        ch1(1,ki,j) = sn*(cc1(1,ki,j)-cc1(2,ki,jc))
        ch1(2,ki,j) = sn*(cc1(2,ki,j)+cc1(1,ki,jc))
        ch1(1,ki,jc) = sn*(cc1(1,ki,j)+cc1(2,ki,jc))
        ch1(2,ki,jc) = sn*(cc1(2,ki,j)-cc1(1,ki,jc))
      end do
    end do

  else

    sn = 1.0E+00 / real ( ip * l1 )

    do ki = 1, lid
      cc1(1,ki,1) = sn * cc1(1,ki,1)
      cc1(2,ki,1) = sn * cc1(2,ki,1)
    end do

    do j = 2, ipph
      jc = ipp2 - j
      do ki = 1, lid
        chold1 = sn*(cc1(1,ki,j)-cc1(2,ki,jc))
        chold2 = sn*(cc1(1,ki,j)+cc1(2,ki,jc))
        cc1(1,ki,j) = chold1
        cc1(2,ki,jc) = sn*(cc1(2,ki,j)-cc1(1,ki,jc))
        cc1(2,ki,j) = sn*(cc1(2,ki,j)+cc1(1,ki,jc))
        cc1(1,ki,jc) = chold2
      end do
    end do

  end if

  return
end
subroutine c1fm1b ( n, inc, c, ch, wa, fnf, fac )

!*****************************************************************************80
!
!! C1FM1B is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  complex c(*)
  real ch(*)
  real fac(*)
  real fnf
  integer ido
  integer inc
  integer inc2
  integer ip
  integer iw
  integer k
  integer k1
  integer l1
  integer l2
  integer lid
  integer n
  integer na
  integer nbr
  integer nf
  real wa(*)

  inc2 = inc + inc
  nf = int ( fnf )
  na = 0
  l1 = 1
  iw = 1

  do k1 = 1, nf

    ip = int ( fac(k1) )
    l2 = ip * l1
    ido = n / l2
    lid = l1 * ido
    nbr = 1 + na + 2 * min ( ip - 2, 4 )

    if ( nbr == 1 ) then
      call c1f2kb ( ido, l1, na, c, inc2, ch, 2, wa(iw) )
    else if ( nbr == 2 ) then
      call c1f2kb ( ido, l1, na, ch, 2, c, inc2, wa(iw) )
    else if ( nbr == 3 ) then
      call c1f3kb ( ido, l1, na, c, inc2, ch, 2, wa(iw) )
    else if ( nbr == 4 ) then
      call c1f3kb ( ido, l1, na, ch, 2, c, inc2, wa(iw) )
    else if ( nbr == 5 ) then
      call c1f4kb ( ido, l1, na, c, inc2, ch, 2, wa(iw) )
    else if ( nbr == 6 ) then
      call c1f4kb ( ido, l1, na, ch, 2, c, inc2, wa(iw) )
    else if ( nbr == 7 ) then
      call c1f5kb ( ido, l1, na, c, inc2, ch, 2, wa(iw) )
    else if ( nbr == 8 ) then
      call c1f5kb ( ido, l1, na, ch, 2, c, inc2, wa(iw) )
    else if ( nbr == 9 ) then
      call c1fgkb ( ido, ip, l1, lid, na, c, c, inc2, ch, ch, 2, wa(iw) )
    else if ( nbr == 10 ) then
      call c1fgkb ( ido, ip, l1, lid, na, ch, ch, 2, c, c, inc2, wa(iw) )
    end if

    l1 = l2
    iw = iw + ( ip - 1 ) * ( ido + ido )

    if ( ip <= 5 ) then
      na = 1 - na
    end if

  end do

  return
end
subroutine c1fm1f ( n, inc, c, ch, wa, fnf, fac )

!*****************************************************************************80
!
!! C1FM1F is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  complex c(*)
  real ch(*)
  real fac(*)
  real fnf
  integer ido
  integer inc
  integer inc2
  integer ip
  integer iw
  integer k1
  integer l1
  integer l2
  integer lid
  integer n
  integer na
  integer nbr
  integer nf
  real wa(*)

  inc2 = inc + inc
  nf = int ( fnf )
  na = 0
  l1 = 1
  iw = 1

  do k1 = 1, nf

     ip = int ( fac(k1) )
     l2 = ip * l1
     ido = n / l2
     lid = l1 * ido
     nbr = 1 + na + 2 * min ( ip - 2, 4 )

     if ( nbr == 1 ) then
       call c1f2kf ( ido, l1, na, c, inc2, ch, 2, wa(iw) )
     else if ( nbr == 2 ) then
       call c1f2kf ( ido, l1, na, ch, 2, c, inc2, wa(iw) )
     else if ( nbr == 3 ) then
       call c1f3kf ( ido, l1, na, c, inc2, ch, 2, wa(iw) )
     else if ( nbr == 4 ) then
       call c1f3kf ( ido, l1, na, ch, 2, c, inc2, wa(iw) )
     else if ( nbr == 5 ) then
       call c1f4kf ( ido, l1, na, c, inc2, ch, 2, wa(iw) )
     else if ( nbr == 6 ) then
       call c1f4kf ( ido, l1, na, ch, 2, c, inc2, wa(iw) )
     else if ( nbr == 7 ) then
       call c1f5kf ( ido, l1, na, c, inc2, ch, 2, wa(iw) )
     else if ( nbr == 8 ) then
       call c1f5kf ( ido, l1, na, ch, 2, c, inc2, wa(iw) )
     else if ( nbr == 9 ) then
       call c1fgkf ( ido, ip, l1, lid, na, c, c, inc2, ch, ch, 1, wa(iw) )
     else if ( nbr == 10 ) then
       call c1fgkf ( ido, ip, l1, lid, na, ch, ch, 2, c, c, inc2, wa(iw) )
     end if

     l1 = l2
     iw = iw + ( ip - 1 ) * ( ido + ido )

     if ( ip <= 5 ) then
       na = 1 - na
     end if

  end do

  return
end
subroutine cfft1b ( n, inc, c, lenc, wsave, lensav, work, lenwrk, ier )

!*****************************************************************************80
!
!! CFFT1B: complex backward fast Fourier transform, 1D.
!
!  Discussion:
!
!    CFFT1B computes the one-dimensional Fourier transform of a single 
!    periodic sequence within a complex array.  This transform is referred 
!    to as the backward transform or Fourier synthesis, transforming the
!    sequence from spectral to physical space.
!
!    This transform is normalized since a call to CFFT1B followed
!    by a call to CFFT1F (or vice-versa) reproduces the original
!    array within roundoff error.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Modified:
!
!    22 March 2005
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
!    Input, integer N, the length of the sequence to be transformed.  The 
!    transform is most efficient when N is a product of small primes.
!
!    Input, integer INC, the increment between the locations, in array C, 
!    of two consecutive elements within the sequence to be transformed.
!
!    Input/output, complex C(LENC) containing the sequence to be transformed.
!
!    Input, integer LENC, the dimension of the C array.  LENC must be at least 
!    INC*(N-1) + 1.
!
!    Input, real WSAVE(LENSAV).  WSAVE's contents must be initialized with 
!    a call to CFFT1I before the first call to routine CFFT1F or CFFT1B 
!    for a given transform length N.  WSAVE's contents may be re-used for 
!    subsequent calls to CFFT1F and CFFT1B with the same N.
!
!    Input, integer LENSAV, the dimension of the WSAVE array.  LENSAV must be
!    at least 2*N + INT(LOG(REAL(N))) + 4.
!
!    Workspace, real WORK(LENWRK).
!
!    Input, integer LENWRK, the dimension of the WORK array.  LENWRK must be at
!    least 2*N.
!
!    Output, integer IER, error flag.
!    0, successful exit;
!    1, input parameter LENC not big enough;
!    2, input parameter LENSAV not big enough;
!    3, input parameter LENWRK not big enough;
!    20, input error returned by lower level routine.
!
  implicit none

  integer lenc
  integer lensav
  integer lenwrk

  complex c(lenc)
  integer ier
  integer inc
  integer iw1
  integer n
  real work(lenwrk)
  real wsave(lensav)

  ier = 0

  if ( lenc < inc * ( n - 1 ) + 1 ) then
    ier = 1
    call xerfft ( 'CFFT1B', 6 )
    return
  end if

  if ( lensav < 2 * n + int ( log ( real ( n ) ) ) + 4 ) then
    ier = 2
    call xerfft ( 'CFFT1B', 8 )
    return
  end if

  if ( lenwrk < 2 * n ) then
    ier = 3
    call xerfft ( 'CFFT1B', 10 )
    return
  end if

  if ( n == 1 ) then
    return
  end if

  iw1 = n + n + 1

  call c1fm1b ( n, inc, c, work, wsave, wsave(iw1), wsave(iw1+1) )

  return
end
subroutine cfft1f ( n, inc, c, lenc, wsave, lensav, work, lenwrk, ier )

!*****************************************************************************80
!
!! CFFT1F: complex forward fast Fourier transform, 1D.
!
!  Discussion:
!
!    CFFT1F computes the one-dimensional Fourier transform of a single 
!    periodic sequence within a complex array.  This transform is referred 
!    to as the forward transform or Fourier analysis, transforming the 
!    sequence from physical to spectral space.
!
!    This transform is normalized since a call to CFFT1F followed
!    by a call to CFFT1B (or vice-versa) reproduces the original
!    array within roundoff error.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Modified:
!
!    22 March 2005
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
!    Input, integer N, the length of the sequence to be transformed.  The 
!    transform is most efficient when N is a product of small primes.
!
!    Input, integer INC, the increment between the locations, in array C, 
!    of two consecutive elements within the sequence to be transformed.
!
!    Input/output, complex C(LENC) containing the sequence to be transformed.
!
!    Input, integer LENC, the dimension of the C array.  LENC must be at least 
!    INC*(N-1) + 1.
!
!    Input, real WSAVE(LENSAV).  WSAVE's contents must be initialized with 
!    a call to CFFT1I before the first call to routine CFFT1F or CFFT1B 
!    for a given transform length N.  WSAVE's contents may be re-used for 
!    subsequent calls to CFFT1F and CFFT1B with the same N.
!
!    Input, integer LENSAV, the dimension of the WSAVE array.  LENSAV must be
!    at least 2*N + INT(LOG(REAL(N))) + 4.
!
!    Workspace, real WORK(LENWRK).
!
!    Input, integer LENWRK, the dimension of the WORK array.  LENWRK must be at
!    least 2*N.
!
!    Output, integer IER, error flag.
!    0, successful exit;
!    1, input parameter LENC   not big enough;
!    2, input parameter LENSAV not big enough;
!    3, input parameter LENWRK not big enough;
!    20, input error returned by lower level routine.
!
  implicit none

  integer lenc
  integer lensav
  integer lenwrk

  complex c(lenc)
  integer ier
  integer inc
  integer iw1
  integer n
  real work(lenwrk)
  real wsave(lensav)

  ier = 0

  if ( lenc < inc * ( n - 1 ) + 1 ) then
    ier = 1
    call xerfft ( 'CFFT1F', 6 )
    return
  end if

  if ( lensav < 2 * n + int ( log ( real ( n ) ) ) + 4 ) then
    ier = 2
    call xerfft ( 'CFFT1F', 8 )
    return
  end if

  if ( lenwrk < 2 * n ) then
    ier = 3
    call xerfft ( 'CFFT1F', 10 )
    return
  end if

  if ( n == 1 ) then
    return
  end if

  iw1 = n + n + 1

  call c1fm1f ( n, inc, c, work, wsave, wsave(iw1), wsave(iw1+1) )

  return
end
subroutine cfft1i ( n, wsave, lensav, ier )

!*****************************************************************************80
!
!! CFFT1I: initialization for CFFT1B and CFFT1F.
!
!  Discussion:
!
!    CFFT1I initializes array WSAVE for use in its companion routines 
!    CFFT1B and CFFT1F.  Routine CFFT1I must be called before the first 
!    call to CFFT1B or CFFT1F, and after whenever the value of integer 
!    N changes.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Modified:
!
!    22 March 2005
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
!    Input, integer N, the length of the sequence to be transformed.  The 
!    transform is most efficient when N is a product of small primes.
!
!    Input, integer LENSAV, the dimension of the WSAVE array.  LENSAV must 
!    be at least 2*N + INT(LOG(REAL(N))) + 4.
!
!    Output, real WSAVE(LENSAV), containing the prime factors of N and 
!    also containing certain trigonometric values which will be used in 
!    routines CFFT1B or CFFT1F.
!
!    Output, integer IER, error flag.
!    0, successful exit;
!    2, input parameter LENSAV not big enough.

  implicit none

  integer lensav

  integer ier
  integer iw1
  integer n
  real wsave(lensav)

  ier = 0

  if ( lensav < 2 * n + int ( log ( real ( n ) ) ) + 4 ) then
    ier = 2
    call xerfft ( 'CFFT1I', 3 )
    return
  end if

  if ( n == 1 ) then
    return
  end if

  iw1 = n + n + 1

  call mcfti1 ( n, wsave, wsave(iw1), wsave(iw1+1) )

  return
end
subroutine cfft2b ( ldim, l, m, c, wsave, lensav, work, lenwrk, ier )

!*****************************************************************************80
!
!! CFFT2B: complex backward fast Fourier transform, 2D.
!
!  Discussion:
!
!    CFFT2B computes the two-dimensional discrete Fourier transform of a 
!    complex periodic array.  This transform is known as the backward 
!    transform or Fourier synthesis, transforming from spectral to 
!    physical space.  Routine CFFT2B is normalized, in that a call to 
!    CFFT2B followed by a call to CFFT2F (or vice-versa) reproduces the 
!    original array within roundoff error. 
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Modified:
!
!    23 March 2005
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
!    Input, integer LDIM, the first dimension of C. 
!
!    Input, integer L, the number of elements to be transformed in the 
!    first dimension of the two-dimensional complex array C.  The value 
!    of L must be less than or equal to that of LDIM.  The transform is
!    most efficient when L is a product of small primes.
!
!    Input, integer M, the number of elements to be transformed in the
!    second dimension of the two-dimensional complex array C.  The transform
!    is most efficient when M is a product of small primes. 
!
!    Input/output, complex C(LDIM,M), on intput, the array of two dimensions
!    containing the (L,M) subarray to be transformed.  On output, the
!    transformed data.
!
!    Input, real WSAVE(LENSAV). WSAVE's contents must be initialized with
!    a call to CFFT2I before the first call to routine CFFT2F
!    or CFFT2B with transform lengths L and M.  WSAVE's contents may be
!    re-used for subsequent calls to CFFT2F and CFFT2B with the same 
!    transform lengths L and M. 
!
!    Input, integer LENSAV, the dimension of the WSAVE array. LENSAV must be
!    at least 2*(L+M) + INT(LOG(REAL(L))) + INT(LOG(REAL(M))) + 8. 
!
!    Workspace, real WORK(LENWRK). 
!
!    Input, integer LENWRK, the dimension of the WORK array.  LENWRK must 
!    be at least 2*L*M. 
!
!    Output, integer IER, the error flag.
!    0, successful exit;
!    2, input parameter LENSAV not big enough;
!    3, input parameter LENWRK not big enough;
!    5, input parameter LDIM < L;
!    20, input error returned by lower level routine.
!
  implicit none

  integer m
  integer ldim
  integer lensav
  integer lenwrk

  complex c(ldim,m)
  integer ier
  integer ier1
  integer iw
  integer l
  real work(lenwrk)
  real wsave(lensav)

  ier = 0

  if ( ldim < l ) then
    ier = 5
    call xerfft ( 'CFFT2B', -2 )
    return
  end if

  if ( lensav < 2 * l + int ( log ( real ( l ) ) ) + &
                2 * m + int ( log ( real ( m ) ) ) + 8 ) then
    ier = 2
    call xerfft ( 'CFFT2B', 6 )
    return
  end if

  if ( lenwrk < 2 * l * m ) then
    ier = 3
    call xerfft ( 'CFFT2B', 8 )
    return
  end if
!
!  Transform the X lines of the C array.
!
  iw = 2 * l + int ( log ( real ( l ) ) * log ( 2.0E+00 ) ) + 3

  call cfftmb ( l, 1, m, ldim, c, (l-1)+ldim*(m-1) +1, &
    wsave(iw), 2*m + int(log(real(m))) + 4, work, 2*l*m, ier1 )

  if ( ier1 /= 0 ) then
    ier = 20
    call xerfft ( 'CFFT2B', -5 )
    return
  end if
!
!  Transform the Y lines of the C array.
!
  iw = 1
  call cfftmb ( m, ldim, l, 1, c, (m-1)*ldim + l, wsave(iw), &
    2*l + int(log(real(l))) + 4, work, 2*m*l, ier1 )

  if ( ier1 /= 0 ) then
    ier = 20
    call xerfft ( 'CFFT2B', -5 )
    return
  end if

  return
end
subroutine cfft2f ( ldim, l, m, c, wsave, lensav, work, lenwrk, ier )

!*****************************************************************************80
!
!! CFFT2F: complex forward fast Fourier transform, 2D.
!
!  Discussion:
!
!    CFFT2F computes the two-dimensional discrete Fourier transform of 
!    a complex periodic array. This transform is known as the forward 
!    transform or Fourier analysis, transforming from physical to 
!    spectral space. Routine CFFT2F is normalized, in that a call to 
!    CFFT2F followed by a call to CFFT2B (or vice-versa) reproduces the 
!    original array within roundoff error. 
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Modified:
!
!    23 March 2005
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
!    Input, integer LDIM, the first dimension of the array C. 
!
!    Input, integer L, the number of elements to be transformed in the first 
!    dimension of the two-dimensional complex array C.  The value of L must 
!    be less than or equal to that of LDIM.  The transform is most efficient 
!    when L is a product of small primes. 
!
!    Input, integer M, the number of elements to be transformed in the 
!    second dimension of the two-dimensional complex array C.  The transform
!    is most efficient when M is a product of small primes. 
!
!    Input/output, complex C(LDIM,M), on input, the array of two dimensions
!    containing the (L,M) subarray to be transformed.  On output, the
!    transformed data.
!
!    Input, real WSAVE(LENSAV). WSAVE's contents must be initialized with a 
!    call to CFFT2I before the first call to routine CFFT2F or 
!    CFFT2B with transform lengths L and M.  WSAVE's contents may be re-used 
!    for subsequent calls to CFFT2F and CFFT2B having those same 
!    transform lengths. 
!
!    Input, integer LENSAV, the dimension of the WSAVE array.  LENSAV must be 
!    at least 2*(L+M) + INT(LOG(REAL(L))) + INT(LOG(REAL(M))) + 8. 
!
!    Workspace, real WORK(LENWRK). 
!
!    Input, integer LENWRK, the dimension of the WORK array.  LENWRK must
!    be at least 2*L*M.
!
!    Output, integer IER, error flag.
!    0, successful exit; 
!    2, input parameter LENSAV not big enough;
!    3, input parameter LENWRK not big enough;
!    5, input parameter LDIM < L;
!    20, input error returned by lower level routine.
!
  implicit none

  integer m
  integer ldim
  integer lensav
  integer lenwrk

  complex c(ldim,m)
  integer ier
  integer ier1
  integer iw
  integer l
  real work(lenwrk)
  real wsave(lensav)

  ier = 0

  if ( ldim < l ) then
    ier = 5
    call xerfft ( 'CFFT2F', -2 )
    return
  end if

  if ( lensav < 2 * l + int ( log ( real ( l ) ) ) + &
                2 * m + int ( log ( real ( m ) ) ) + 8 ) then
    ier = 2
    call xerfft ( 'CFFT2F', 6 )
    return
  end if

  if ( lenwrk < 2 * l * m ) then
    ier = 3
    call xerfft ( 'CFFT2F', 8 )
    return
  end if
!
!  Transform the X lines of the C array.
!
  iw = 2 * l + int ( log ( real ( l ) ) * log ( 2.0E+00 ) ) + 3

  call cfftmf ( l, 1, m, ldim, c, (l-1) + ldim*(m-1) +1, wsave(iw), &
    2*m + int(log(real(m))) + 4, work, 2*l*m, ier1 )

  if ( ier1 /= 0 ) then
    ier = 20
    call xerfft ( 'CFFT2F', -5 )
    return
  end if
!
!  Transform the Y lines of the C array.
!
  iw = 1

  call cfftmf ( m, ldim, l, 1, c, (m-1)*ldim + l, wsave(iw), &
    2*l + int(log(real(l))) + 4, work, 2*m*l, ier1 )

  if ( ier1 /= 0 ) then
    ier = 20
    call xerfft ( 'CFFT2F', -5 )
    return
  end if

  return
end
subroutine cfft2i ( l, m, wsave, lensav, ier )

!*****************************************************************************80
!
!! CFFT2I: initialization for CFFT2B and CFFT2F.
!
!  Discussion:
!
!    CFFT2I initializes real array WSAVE for use in its companion 
!    routines CFFT2F and CFFT2B for computing two-dimensional fast 
!    Fourier transforms of complex data.  Prime factorizations of L and M,
!    together with tabulations of the trigonometric functions, are 
!    computed and stored in array WSAVE.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Modified:
!
!    23 March 2005
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
!    Input, integer L, the number of elements to be transformed in the 
!    first dimension.  The transform is most efficient when L is a product
!    of small primes.
!
!    Input, integer M, the number of elements to be transformed in the 
!    second dimension.  The transform is most efficient when M is a product
!    of small primes. 
!
!    Input, integer LENSAV, the dimension of the WSAVE array.  LENSAV must 
!    be at least 2*(L+M) + INT(LOG(REAL(L))) + INT(LOG(REAL(M))) + 8.
!
!    Output, real WSAVE(LENSAV), contains the prime factors of L and M, and 
!    also certain trigonometric values which will be used in routines 
!    CFFT2B or CFFT2F. 
!
!    Output, integer IER, error flag.
!    0, successful exit;
!    2, input parameter LENSAV not big enough;
!    20, input error returned by lower level routine.
!
  implicit none

  integer lensav

  integer ier
  integer ier1
  integer l
  integer m
  real wsave(lensav)

  ier = 0

  if ( lensav < 2 * l + int ( log ( real ( l ) ) ) + &
                2 * m + int ( log ( real ( m ) ) ) + 8 ) then
    ier = 2
    call xerfft ( 'CFFT2I', 4 )
    return
  end if

  call cfftmi ( l, wsave(1), 2*l + int(log(real(l))) + 4, ier1 )

  if ( ier1 /= 0 ) then
    ier = 20
    call xerfft ( 'CFFT2I', -5 )
    return
  end if

  call cfftmi ( m, wsave(2*l+int(log(real(l))*log( 2.0E+00 )) + 3), &
    2*m + int(log(real(m))) + 4, ier1 )

  if ( ier1 /= 0 ) then
    ier = 20
    call xerfft ( 'CFFT2I', -5 )
    return
  end if

  return
end
subroutine cfftmb ( lot, jump, n, inc, c, lenc, wsave, lensav, work, &
  lenwrk, ier )

!*****************************************************************************80
!
!! CFFTMB: complex backward fast Fourier transform, 1D, multiple vectors.
!
!  Discussion:
!
!    CFFTMB computes the one-dimensional Fourier transform of multiple 
!    periodic sequences within a complex array.  This transform is referred
!    to as the backward transform or Fourier synthesis, transforming the
!    sequences from spectral to physical space.  This transform is 
!    normalized since a call to CFFTMF followed by a call to CFFTMB (or
!    vice-versa) reproduces the original array within roundoff error. 
!
!    The parameters INC, JUMP, N and LOT are consistent if equality 
!    I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT 
!    implies I1=I2 and J1=J2.  For multiple FFTs to execute correctly, 
!    input variables INC, JUMP, N and LOT must be consistent, otherwise 
!    at least one array element mistakenly is transformed more than once.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Modified:
!
!    24 March 2005
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
!    Input, integer LOT, the number of sequences to be transformed within
!    array C. 
!
!    Input, integer JUMP, the increment between the locations, in array C,
!    of the first elements of two consecutive sequences to be transformed. 
!
!    Input, integer N, the length of each sequence to be transformed.  The
!    transform is most efficient when N is a product of small primes.
!
!    Input, integer INC, the increment between the locations, in array C,
!    of two consecutive elements within the same sequence to be transformed. 
!
!    Input/output, complex C(LENC), an array containing LOT sequences,
!    each having length N, to be transformed.  C can have any number of
!    dimensions, but the total number of locations must be at least LENC. 
!    On output, C contains the transformed sequences.
!
!    Input, integer LENC, the dimension of the C array.  LENC must be at least
!    (LOT-1)*JUMP + INC*(N-1) + 1. 
!
!    Input, real WSAVE(LENSAV).  WSAVE's contents must be initialized with 
!    a call to CFFTMI before the first call to routine CFFTMF 
!    or CFFTMB for a given transform length N. 
!
!    Input, integer LENSAV, the dimension of the WSAVE array.  LENSAV must be
!    at least 2*N + INT(LOG(REAL(N))) + 4.
!
!    Workspace, real WORK(LENWRK). 
!
!    Input, integer LENWRK, the dimension of the WORK array.  LENWRK must be at
!    least 2*LOT*N.
!
!    Output, integer IER, error flag.
!    0, successful exit
!    1, input parameter LENC not big enough;
!    2, input parameter LENSAV not big enough;
!    3, input parameter LENWRK not big enough;
!    4, input parameters INC, JUMP, N, LOT are not consistent. 
!
  implicit none

  integer lenc
  integer lensav
  integer lenwrk

  complex c(lenc)
  integer ier
  integer inc
  integer iw1
  integer jump
  integer lot
  integer n
  real work(lenwrk)
  real wsave(lensav)
  logical xercon

  ier = 0

  if ( lenc < ( lot - 1 ) * jump + inc * ( n - 1 ) + 1 ) then
    ier = 1
    call xerfft ( 'CFFTMB', 6 )
    return
  end if

  if ( lensav < 2 * n + int ( log ( real ( n ) ) ) + 4 ) then
    ier = 2
    call xerfft ( 'CFFTMB', 8 )
    return
  end if

  if ( lenwrk < 2 * lot * n ) then
    ier = 3
    call xerfft ( 'CFFTMB', 10 )
    return
  end if

  if ( .not. xercon ( inc, jump, n, lot ) ) then
    ier = 4
    call xerfft ( 'CFFTMB', -1 )
    return
  end if

  if ( n == 1 ) then
    return
  end if

  iw1 = n + n + 1

  call cmfm1b ( lot, jump, n, inc, c, work, wsave, wsave(iw1), &
    wsave(iw1+1) )

  return
end
subroutine cfftmf ( lot, jump, n, inc, c, lenc, wsave, lensav, work, &
  lenwrk, ier )

!*****************************************************************************80
!
!! CFFTMF: complex forward fast Fourier transform, 1D, multiple vectors.
!
!  Discussion:
!
!    CFFTMF computes the one-dimensional Fourier transform of multiple 
!    periodic sequences within a complex array. This transform is referred 
!    to as the forward transform or Fourier analysis, transforming the 
!    sequences from physical to spectral space. This transform is 
!    normalized since a call to CFFTMF followed by a call to CFFTMB 
!    (or vice-versa) reproduces the original array within roundoff error. 
!
!    The parameters integers INC, JUMP, N and LOT are consistent if equality
!    I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT 
!    implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, 
!    input variables INC, JUMP, N and LOT must be consistent, otherwise 
!    at least one array element mistakenly is transformed more than once.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Modified:
!
!    24 March 2005
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
!    Input, integer LOT, the number of sequences to be transformed within
!    array C. 
!
!    Input, integer JUMP, the increment between the locations, in array C,
!    of the first elements of two consecutive sequences to be transformed.
!
!    Input, integer N, the length of each sequence to be transformed. 
!    The transform is most efficient when N is a product of small primes. 
!
!    Input, integer INC, the increment between the locations, in array C, 
!    of two consecutive elements within the same sequence to be transformed.
!
!    Input/output, complex C(LENC), array containing LOT sequences, each 
!    having length N, to be transformed.  C can have any number of dimensions,
!    but the total number of locations must be at least LENC.
!
!    Input, integer LENC, the dimension of the C array.  LENC must be at least
!    (LOT-1)*JUMP + INC*(N-1) + 1. 
!
!    Input, real WSAVE(LENSAV).  WSAVE's contents must be initialized with 
!    a call to CFFTMI before the first call to routine CFFTMF 
!    or CFFTMB for a given transform length N. 
!
!    Input, integer LENSAV, the dimension of the WSAVE array.  LENSAV must be 
!    at least 2*N + INT(LOG(REAL(N))) + 4. 
!
!    Workspace, real WORK(LENWRK). 
!
!    Input, integer LENWRK, the dimension of the WORK array.  LENWRK must be
!    at least 2*LOT*N. 
!
!    Output, integer IER, error flag.
!    0 successful exit;
!    1 input parameter LENC not big enough;
!    2 input parameter LENSAV not big enough;
!    3 input parameter LENWRK not big enough;
!    4 input parameters INC, JUMP, N, LOT are not consistent.
!
  implicit none

  integer lenc
  integer lensav
  integer lenwrk

  complex c(lenc)
  integer ier
  integer inc
  integer iw1
  integer jump
  integer lot
  integer n
  real work(lenwrk)
  real wsave(lensav)
  logical xercon

  ier = 0

  if ( lenc < ( lot - 1 ) * jump + inc * ( n - 1 ) + 1 ) then
    ier = 1
    call xerfft ( 'CFFTMF', 6 )
    return
  end if

  if ( lensav < 2 * n + int ( log ( real ( n ) ) ) + 4 ) then
    ier = 2
    call xerfft ( 'CFFTMF', 8 )
    return
  end if

  if ( lenwrk < 2 * lot * n ) then
    ier = 3
    call xerfft ( 'CFFTMF', 10 )
    return
  end if

  if ( .not. xercon ( inc, jump, n, lot ) ) then
    ier = 4
    call xerfft ( 'CFFTMF', -1 )
    return
  end if

  if ( n == 1 ) then
    return
  end if

  iw1 = n + n + 1

  call cmfm1f ( lot, jump, n, inc, c, work, wsave, wsave(iw1), &
    wsave(iw1+1) )

  return
end
subroutine cfftmi ( n, wsave, lensav, ier )

!*****************************************************************************80
!
!! CFFTMI: initialization for CFFTMB and CFFTMF.
!
!  Discussion:
!
!    CFFTMI initializes array WSAVE for use in its companion routines 
!    CFFTMB and CFFTMF.  CFFTMI must be called before the first call 
!    to CFFTMB or CFFTMF, and after whenever the value of integer N changes. 
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Modified:
!
!    24 March 2005
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
!    Input, integer N, the length of each sequence to be transformed. 
!    The transform is most efficient when N is a product of small primes. 
!
!    Input, integer LENSAV, the dimension of the WSAVE array.  LENSAV must be
!    at least 2*N + INT(LOG(REAL(N))) + 4. 
!
!    Output, real WSAVE(LENSAV), containing the prime factors of N and
!    also containing certain trigonometric values which will be used in
!    routines CFFTMB or CFFTMF.
!
!    Output, integer IER, error flag.
!    0, successful exit; 
!    2, input parameter LENSAV not big enough.
!
  implicit none

  integer lensav

  integer ier
  integer iw1
  integer n
  real wsave(lensav)

  ier = 0

  if ( lensav < 2 * n + int ( log ( real ( n ) ) ) + 4 ) then
    ier = 2
    call xerfft ( 'cfftmi ', 3 )
    return
  end if

  if ( n == 1 ) then
    return
  end if

  iw1 = n + n + 1
  call mcfti1 ( n, wsave, wsave(iw1), wsave(iw1+1) )

  return
end
subroutine cmf2kb ( lot, ido, l1, na, cc, im1, in1, ch, im2, in2, wa )

!*****************************************************************************80
!
!! CMF2KB is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(2,in1,l1,ido,2)
  real ch(2,in2,l1,2,ido)
  real chold1
  real chold2
  integer i
  integer im1
  integer im2
  integer k
  integer lot
  integer m1
  integer m1d
  integer m2
  integer m2s
  integer na
  real ti2
  real tr2
  real wa(ido,1,2)

  m1d = ( lot - 1 ) * im1 + 1
  m2s = 1 - im2

  if ( 1 < ido .or. na == 1 ) then

    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch(1,m2,k,1,1) = cc(1,m1,k,1,1) + cc(1,m1,k,1,2)
        ch(1,m2,k,2,1) = cc(1,m1,k,1,1) - cc(1,m1,k,1,2)
        ch(2,m2,k,1,1) = cc(2,m1,k,1,1) + cc(2,m1,k,1,2)
        ch(2,m2,k,2,1) = cc(2,m1,k,1,1) - cc(2,m1,k,1,2)
      end do
    end do

    do i = 2, ido
      do k = 1, l1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ch(1,m2,k,1,i) = cc(1,m1,k,i,1)+cc(1,m1,k,i,2)
          tr2 = cc(1,m1,k,i,1)-cc(1,m1,k,i,2)
          ch(2,m2,k,1,i) = cc(2,m1,k,i,1)+cc(2,m1,k,i,2)
          ti2 = cc(2,m1,k,i,1)-cc(2,m1,k,i,2)

          ch(2,m2,k,2,i) = wa(i,1,1) * ti2 + wa(i,1,2) * tr2
          ch(1,m2,k,2,i) = wa(i,1,1) * tr2 - wa(i,1,2) * ti2

        end do
      end do
    end do

  else

    do k = 1, l1
      do m1 = 1, m1d, im1

        chold1         = cc(1,m1,k,1,1) + cc(1,m1,k,1,2)
        cc(1,m1,k,1,2) = cc(1,m1,k,1,1) - cc(1,m1,k,1,2)
        cc(1,m1,k,1,1) = chold1

        chold2         = cc(2,m1,k,1,1) + cc(2,m1,k,1,2)
        cc(2,m1,k,1,2) = cc(2,m1,k,1,1) - cc(2,m1,k,1,2)
        cc(2,m1,k,1,1) = chold2

      end do
    end do

  end if

  return
end
subroutine cmf2kf ( lot, ido, l1, na, cc, im1, in1, ch, im2, in2, wa )

!*****************************************************************************80
!
!! CMF2KF is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(2,in1,l1,ido,2)
  real ch(2,in2,l1,2,ido)
  real chold1
  real chold2
  integer i
  integer im1
  integer im2
  integer k
  integer lid
  integer lot
  integer m1
  integer m1d
  integer m2
  integer m2s
  integer n
  integer na
  real sn
  real ti2
  real tr2
  real wa(ido,1,2)

  m1d = ( lot - 1 ) * im1 + 1
  m2s = 1 - im2

  if ( 1 < ido ) then

    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch(1,m2,k,1,1) = cc(1,m1,k,1,1)+cc(1,m1,k,1,2)
        ch(1,m2,k,2,1) = cc(1,m1,k,1,1)-cc(1,m1,k,1,2)
        ch(2,m2,k,1,1) = cc(2,m1,k,1,1)+cc(2,m1,k,1,2)
        ch(2,m2,k,2,1) = cc(2,m1,k,1,1)-cc(2,m1,k,1,2)
      end do
    end do

    do i = 2, ido
      do k = 1, l1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ch(1,m2,k,1,i) = cc(1,m1,k,i,1)+cc(1,m1,k,i,2)
          tr2 = cc(1,m1,k,i,1)-cc(1,m1,k,i,2)
          ch(2,m2,k,1,i) = cc(2,m1,k,i,1)+cc(2,m1,k,i,2)
          ti2 = cc(2,m1,k,i,1)-cc(2,m1,k,i,2)
          ch(2,m2,k,2,i) = wa(i,1,1)*ti2-wa(i,1,2)*tr2
          ch(1,m2,k,2,i) = wa(i,1,1)*tr2+wa(i,1,2)*ti2
        end do
      end do
    end do

  else if ( na == 1 ) then

    sn = 1.0E+00 / real ( 2 * l1 )

    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch(1,m2,k,1,1) = sn * ( cc(1,m1,k,1,1) + cc(1,m1,k,1,2) )
        ch(1,m2,k,2,1) = sn * ( cc(1,m1,k,1,1) - cc(1,m1,k,1,2) )
        ch(2,m2,k,1,1) = sn * ( cc(2,m1,k,1,1) + cc(2,m1,k,1,2) )
        ch(2,m2,k,2,1) = sn * ( cc(2,m1,k,1,1) - cc(2,m1,k,1,2) )
      end do
    end do

  else

    sn = 1.0E+00 / real ( 2 * l1 )

    do k = 1, l1
      do m1 = 1, m1d, im1

        chold1         = sn * ( cc(1,m1,k,1,1) + cc(1,m1,k,1,2) )
        cc(1,m1,k,1,2) = sn * ( cc(1,m1,k,1,1) - cc(1,m1,k,1,2) )
        cc(1,m1,k,1,1) = chold1

        chold2         = sn * ( cc(2,m1,k,1,1) + cc(2,m1,k,1,2) )
        cc(2,m1,k,1,2) = sn * ( cc(2,m1,k,1,1) - cc(2,m1,k,1,2) )
        cc(2,m1,k,1,1) = chold2

      end do
    end do

  end if

  return
end
subroutine cmf3kb ( lot, ido, l1, na, cc, im1, in1, ch, im2, in2, wa )

!*****************************************************************************80
!
!! CMF3KB is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(2,in1,l1,ido,3)
  real ch(2,in2,l1,3,ido)
  real ci2
  real ci3
  real cr2
  real cr3
  real di2
  real di3
  real dr2
  real dr3
  integer i
  integer im1
  integer im2
  integer k
  integer lot
  integer m1
  integer m1d
  integer m2
  integer m2s
  integer na
  real, parameter :: taui =  0.866025403784439E+00
  real, parameter :: taur = -0.5E+00
  real ti2
  real tr2
  real wa(ido,2,2)

  m1d = ( lot - 1 ) * im1 + 1
  m2s = 1 - im2

  if ( 1 < ido .or. na == 1 ) then
    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1

        m2 = m2 + im2

        tr2 = cc(1,m1,k,1,2)+cc(1,m1,k,1,3)
        cr2 = cc(1,m1,k,1,1)+taur*tr2
        ch(1,m2,k,1,1) = cc(1,m1,k,1,1)+tr2

        ti2 = cc(2,m1,k,1,2)+cc(2,m1,k,1,3)
        ci2 = cc(2,m1,k,1,1)+taur*ti2
        ch(2,m2,k,1,1) = cc(2,m1,k,1,1)+ti2

        cr3 = taui * (cc(1,m1,k,1,2)-cc(1,m1,k,1,3))
        ci3 = taui * (cc(2,m1,k,1,2)-cc(2,m1,k,1,3))

        ch(1,m2,k,2,1) = cr2-ci3
        ch(1,m2,k,3,1) = cr2+ci3
        ch(2,m2,k,2,1) = ci2+cr3
        ch(2,m2,k,3,1) = ci2-cr3

      end do
    end do

    do i = 2, ido
      do k = 1, l1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          tr2 = cc(1,m1,k,i,2)+cc(1,m1,k,i,3)
          cr2 = cc(1,m1,k,i,1)+taur*tr2
          ch(1,m2,k,1,i) = cc(1,m1,k,i,1)+tr2
          ti2 = cc(2,m1,k,i,2)+cc(2,m1,k,i,3)
          ci2 = cc(2,m1,k,i,1)+taur*ti2
          ch(2,m2,k,1,i) = cc(2,m1,k,i,1)+ti2
          cr3 = taui*(cc(1,m1,k,i,2)-cc(1,m1,k,i,3))
          ci3 = taui*(cc(2,m1,k,i,2)-cc(2,m1,k,i,3))
          dr2 = cr2-ci3
          dr3 = cr2+ci3
          di2 = ci2+cr3
          di3 = ci2-cr3
          ch(2,m2,k,2,i) = wa(i,1,1)*di2+wa(i,1,2)*dr2
          ch(1,m2,k,2,i) = wa(i,1,1)*dr2-wa(i,1,2)*di2
          ch(2,m2,k,3,i) = wa(i,2,1)*di3+wa(i,2,2)*dr3
          ch(1,m2,k,3,i) = wa(i,2,1)*dr3-wa(i,2,2)*di3
        end do
      end do
    end do

  else

    do k = 1, l1
      do m1 = 1, m1d, im1
        tr2 = cc(1,m1,k,1,2)+cc(1,m1,k,1,3)
        cr2 = cc(1,m1,k,1,1)+taur*tr2
        cc(1,m1,k,1,1) = cc(1,m1,k,1,1)+tr2
        ti2 = cc(2,m1,k,1,2)+cc(2,m1,k,1,3)
        ci2 = cc(2,m1,k,1,1)+taur*ti2
        cc(2,m1,k,1,1) = cc(2,m1,k,1,1)+ti2
        cr3 = taui*(cc(1,m1,k,1,2)-cc(1,m1,k,1,3))
        ci3 = taui*(cc(2,m1,k,1,2)-cc(2,m1,k,1,3))
        cc(1,m1,k,1,2) = cr2-ci3
        cc(1,m1,k,1,3) = cr2+ci3
        cc(2,m1,k,1,2) = ci2+cr3
        cc(2,m1,k,1,3) = ci2-cr3
      end do
    end do

  end if

  return
end
subroutine cmf3kf ( lot, ido, l1, na, cc, im1, in1, ch, im2, in2, wa )

!*****************************************************************************80
!
!! CMF3KF is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(2,in1,l1,ido,3)
  real ch(2,in2,l1,3,ido)
  real ci2
  real ci3
  real cr2
  real cr3
  real di2
  real di3
  real dr2
  real dr3
  integer i
  integer im1
  integer im2
  integer k
  integer lot
  integer m1
  integer m1d
  integer m2
  integer m2s
  integer na
  real sn
  real, parameter :: taui = -0.866025403784439E+00
  real, parameter :: taur = -0.5E+00
  real ti2
  real tr2
  real wa(ido,2,2)

  m1d = ( lot - 1 ) * im1 + 1
  m2s = 1 - im2

  if ( 1 < ido ) then

    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        tr2 = cc(1,m1,k,1,2)+cc(1,m1,k,1,3)
        cr2 = cc(1,m1,k,1,1)+taur*tr2
        ch(1,m2,k,1,1) = cc(1,m1,k,1,1)+tr2
        ti2 = cc(2,m1,k,1,2)+cc(2,m1,k,1,3)
        ci2 = cc(2,m1,k,1,1)+taur*ti2
        ch(2,m2,k,1,1) = cc(2,m1,k,1,1)+ti2
        cr3 = taui*(cc(1,m1,k,1,2)-cc(1,m1,k,1,3))
        ci3 = taui*(cc(2,m1,k,1,2)-cc(2,m1,k,1,3))
        ch(1,m2,k,2,1) = cr2-ci3
        ch(1,m2,k,3,1) = cr2+ci3
        ch(2,m2,k,2,1) = ci2+cr3
        ch(2,m2,k,3,1) = ci2-cr3
      end do
    end do

    do i = 2, ido
      do k = 1, l1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          tr2 = cc(1,m1,k,i,2)+cc(1,m1,k,i,3)
          cr2 = cc(1,m1,k,i,1)+taur*tr2
          ch(1,m2,k,1,i) = cc(1,m1,k,i,1)+tr2
          ti2 = cc(2,m1,k,i,2)+cc(2,m1,k,i,3)
          ci2 = cc(2,m1,k,i,1)+taur*ti2
          ch(2,m2,k,1,i) = cc(2,m1,k,i,1)+ti2
          cr3 = taui*(cc(1,m1,k,i,2)-cc(1,m1,k,i,3))
          ci3 = taui*(cc(2,m1,k,i,2)-cc(2,m1,k,i,3))
          dr2 = cr2-ci3
          dr3 = cr2+ci3
          di2 = ci2+cr3
          di3 = ci2-cr3
          ch(2,m2,k,2,i) = wa(i,1,1)*di2-wa(i,1,2)*dr2
          ch(1,m2,k,2,i) = wa(i,1,1)*dr2+wa(i,1,2)*di2
          ch(2,m2,k,3,i) = wa(i,2,1)*di3-wa(i,2,2)*dr3
          ch(1,m2,k,3,i) = wa(i,2,1)*dr3+wa(i,2,2)*di3
        end do
      end do
    end do

  else if ( na == 1 ) then

    sn = 1.0E+00 / real ( 3 * l1 )

    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        tr2 = cc(1,m1,k,1,2)+cc(1,m1,k,1,3)
        cr2 = cc(1,m1,k,1,1)+taur*tr2
        ch(1,m2,k,1,1) = sn*(cc(1,m1,k,1,1)+tr2)
        ti2 = cc(2,m1,k,1,2)+cc(2,m1,k,1,3)
        ci2 = cc(2,m1,k,1,1)+taur*ti2
        ch(2,m2,k,1,1) = sn*(cc(2,m1,k,1,1)+ti2)
        cr3 = taui*(cc(1,m1,k,1,2)-cc(1,m1,k,1,3))
        ci3 = taui*(cc(2,m1,k,1,2)-cc(2,m1,k,1,3))
        ch(1,m2,k,2,1) = sn*(cr2-ci3)
        ch(1,m2,k,3,1) = sn*(cr2+ci3)
        ch(2,m2,k,2,1) = sn*(ci2+cr3)
        ch(2,m2,k,3,1) = sn*(ci2-cr3)
      end do
    end do

  else

    sn = 1.0E+00 / real ( 3 * l1 )

    do k = 1, l1
      do m1 = 1, m1d, im1
        tr2 = cc(1,m1,k,1,2)+cc(1,m1,k,1,3)
        cr2 = cc(1,m1,k,1,1)+taur*tr2
        cc(1,m1,k,1,1) = sn*(cc(1,m1,k,1,1)+tr2)
        ti2 = cc(2,m1,k,1,2)+cc(2,m1,k,1,3)
        ci2 = cc(2,m1,k,1,1)+taur*ti2
        cc(2,m1,k,1,1) = sn*(cc(2,m1,k,1,1)+ti2)
        cr3 = taui*(cc(1,m1,k,1,2)-cc(1,m1,k,1,3))
        ci3 = taui*(cc(2,m1,k,1,2)-cc(2,m1,k,1,3))
        cc(1,m1,k,1,2) = sn*(cr2-ci3)
        cc(1,m1,k,1,3) = sn*(cr2+ci3)
        cc(2,m1,k,1,2) = sn*(ci2+cr3)
        cc(2,m1,k,1,3) = sn*(ci2-cr3)
      end do
    end do

  end if

  return
end
subroutine cmf4kb ( lot, ido, l1, na, cc, im1, in1, ch, im2, in2, wa )

!*****************************************************************************80
!
!! CMF4KB is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(2,in1,l1,ido,4)
  real ch(2,in2,l1,4,ido)
  real ci2
  real ci3
  real ci4
  real cr2
  real cr3
  real cr4
  integer i
  integer im1
  integer im2
  integer k
  integer lot
  integer m1
  integer m1d
  integer m2
  integer m2s
  integer na
  real ti1
  real ti2
  real ti3
  real ti4
  real tr1
  real tr2
  real tr3
  real tr4
  real wa(ido,3,2)

  m1d = ( lot - 1 ) * im1 + 1
  m2s = 1 - im2

  if ( 1 < ido .or. na == 1 ) then

    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ti1 = cc(2,m1,k,1,1)-cc(2,m1,k,1,3)
        ti2 = cc(2,m1,k,1,1)+cc(2,m1,k,1,3)
        tr4 = cc(2,m1,k,1,4)-cc(2,m1,k,1,2)
        ti3 = cc(2,m1,k,1,2)+cc(2,m1,k,1,4)
        tr1 = cc(1,m1,k,1,1)-cc(1,m1,k,1,3)
        tr2 = cc(1,m1,k,1,1)+cc(1,m1,k,1,3)
        ti4 = cc(1,m1,k,1,2)-cc(1,m1,k,1,4)
        tr3 = cc(1,m1,k,1,2)+cc(1,m1,k,1,4)
        ch(1,m2,k,1,1) = tr2+tr3
        ch(1,m2,k,3,1) = tr2-tr3
        ch(2,m2,k,1,1) = ti2+ti3
        ch(2,m2,k,3,1) = ti2-ti3
        ch(1,m2,k,2,1) = tr1+tr4
        ch(1,m2,k,4,1) = tr1-tr4
        ch(2,m2,k,2,1) = ti1+ti4
        ch(2,m2,k,4,1) = ti1-ti4
      end do
    end do

    do i = 2, ido
      do k = 1, l1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ti1 = cc(2,m1,k,i,1)-cc(2,m1,k,i,3)
          ti2 = cc(2,m1,k,i,1)+cc(2,m1,k,i,3)
          ti3 = cc(2,m1,k,i,2)+cc(2,m1,k,i,4)
          tr4 = cc(2,m1,k,i,4)-cc(2,m1,k,i,2)
          tr1 = cc(1,m1,k,i,1)-cc(1,m1,k,i,3)
          tr2 = cc(1,m1,k,i,1)+cc(1,m1,k,i,3)
          ti4 = cc(1,m1,k,i,2)-cc(1,m1,k,i,4)
          tr3 = cc(1,m1,k,i,2)+cc(1,m1,k,i,4)
          ch(1,m2,k,1,i) = tr2+tr3
          cr3 = tr2-tr3
          ch(2,m2,k,1,i) = ti2+ti3
          ci3 = ti2-ti3
          cr2 = tr1+tr4
          cr4 = tr1-tr4
          ci2 = ti1+ti4
          ci4 = ti1-ti4
          ch(1,m2,k,2,i) = wa(i,1,1)*cr2-wa(i,1,2)*ci2
          ch(2,m2,k,2,i) = wa(i,1,1)*ci2+wa(i,1,2)*cr2
          ch(1,m2,k,3,i) = wa(i,2,1)*cr3-wa(i,2,2)*ci3
          ch(2,m2,k,3,i) = wa(i,2,1)*ci3+wa(i,2,2)*cr3
          ch(1,m2,k,4,i) = wa(i,3,1)*cr4-wa(i,3,2)*ci4
          ch(2,m2,k,4,i) = wa(i,3,1)*ci4+wa(i,3,2)*cr4
        end do
      end do
    end do

  else

    do k = 1, l1
      do m1 = 1, m1d, im1
        ti1 = cc(2,m1,k,1,1)-cc(2,m1,k,1,3)
        ti2 = cc(2,m1,k,1,1)+cc(2,m1,k,1,3)
        tr4 = cc(2,m1,k,1,4)-cc(2,m1,k,1,2)
        ti3 = cc(2,m1,k,1,2)+cc(2,m1,k,1,4)
        tr1 = cc(1,m1,k,1,1)-cc(1,m1,k,1,3)
        tr2 = cc(1,m1,k,1,1)+cc(1,m1,k,1,3)
        ti4 = cc(1,m1,k,1,2)-cc(1,m1,k,1,4)
        tr3 = cc(1,m1,k,1,2)+cc(1,m1,k,1,4)
        cc(1,m1,k,1,1) = tr2+tr3
        cc(1,m1,k,1,3) = tr2-tr3
        cc(2,m1,k,1,1) = ti2+ti3
        cc(2,m1,k,1,3) = ti2-ti3
        cc(1,m1,k,1,2) = tr1+tr4
        cc(1,m1,k,1,4) = tr1-tr4
        cc(2,m1,k,1,2) = ti1+ti4
        cc(2,m1,k,1,4) = ti1-ti4
      end do
    end do

  end if

  return
end
subroutine cmf4kf ( lot, ido, l1, na, cc, im1, in1, ch, im2, in2, wa )

!*****************************************************************************80
!
!! CMF4KF is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(2,in1,l1,ido,4)
  real ch(2,in2,l1,4,ido)
  real ci2
  real ci3
  real ci4
  real cr2
  real cr3
  real cr4
  integer i
  integer im1
  integer im2
  integer k
  integer lot
  integer m1
  integer m1d
  integer m2
  integer m2s
  integer na
  real sn
  real ti1
  real ti2
  real ti3
  real ti4
  real tr1
  real tr2
  real tr3
  real tr4
  real wa(ido,3,2)

  m1d = ( lot - 1 ) * im1 + 1
  m2s = 1 - im2

  if ( 1 < ido ) then

    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ti1 = cc(2,m1,k,1,1)-cc(2,m1,k,1,3)
        ti2 = cc(2,m1,k,1,1)+cc(2,m1,k,1,3)
        tr4 = cc(2,m1,k,1,2)-cc(2,m1,k,1,4)
        ti3 = cc(2,m1,k,1,2)+cc(2,m1,k,1,4)
        tr1 = cc(1,m1,k,1,1)-cc(1,m1,k,1,3)
        tr2 = cc(1,m1,k,1,1)+cc(1,m1,k,1,3)
        ti4 = cc(1,m1,k,1,4)-cc(1,m1,k,1,2)
        tr3 = cc(1,m1,k,1,2)+cc(1,m1,k,1,4)
        ch(1,m2,k,1,1) = tr2+tr3
        ch(1,m2,k,3,1) = tr2-tr3
        ch(2,m2,k,1,1) = ti2+ti3
        ch(2,m2,k,3,1) = ti2-ti3
        ch(1,m2,k,2,1) = tr1+tr4
        ch(1,m2,k,4,1) = tr1-tr4
        ch(2,m2,k,2,1) = ti1+ti4
        ch(2,m2,k,4,1) = ti1-ti4
      end do
    end do

    do i = 2, ido
      do k = 1, l1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ti1 = cc(2,m1,k,i,1)-cc(2,m1,k,i,3)
          ti2 = cc(2,m1,k,i,1)+cc(2,m1,k,i,3)
          ti3 = cc(2,m1,k,i,2)+cc(2,m1,k,i,4)
          tr4 = cc(2,m1,k,i,2)-cc(2,m1,k,i,4)
          tr1 = cc(1,m1,k,i,1)-cc(1,m1,k,i,3)
          tr2 = cc(1,m1,k,i,1)+cc(1,m1,k,i,3)
          ti4 = cc(1,m1,k,i,4)-cc(1,m1,k,i,2)
          tr3 = cc(1,m1,k,i,2)+cc(1,m1,k,i,4)
          ch(1,m2,k,1,i) = tr2+tr3
          cr3 = tr2-tr3
          ch(2,m2,k,1,i) = ti2+ti3
          ci3 = ti2-ti3
          cr2 = tr1+tr4
          cr4 = tr1-tr4
          ci2 = ti1+ti4
          ci4 = ti1-ti4
          ch(1,m2,k,2,i) = wa(i,1,1)*cr2+wa(i,1,2)*ci2
          ch(2,m2,k,2,i) = wa(i,1,1)*ci2-wa(i,1,2)*cr2
          ch(1,m2,k,3,i) = wa(i,2,1)*cr3+wa(i,2,2)*ci3
          ch(2,m2,k,3,i) = wa(i,2,1)*ci3-wa(i,2,2)*cr3
          ch(1,m2,k,4,i) = wa(i,3,1)*cr4+wa(i,3,2)*ci4
          ch(2,m2,k,4,i) = wa(i,3,1)*ci4-wa(i,3,2)*cr4
        end do
      end do
    end do

  else if ( na == 1 ) then

    sn = 1.0E+00 / real ( 4 * l1 )

    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ti1 = cc(2,m1,k,1,1)-cc(2,m1,k,1,3)
        ti2 = cc(2,m1,k,1,1)+cc(2,m1,k,1,3)
        tr4 = cc(2,m1,k,1,2)-cc(2,m1,k,1,4)
        ti3 = cc(2,m1,k,1,2)+cc(2,m1,k,1,4)
        tr1 = cc(1,m1,k,1,1)-cc(1,m1,k,1,3)
        tr2 = cc(1,m1,k,1,1)+cc(1,m1,k,1,3)
        ti4 = cc(1,m1,k,1,4)-cc(1,m1,k,1,2)
        tr3 = cc(1,m1,k,1,2)+cc(1,m1,k,1,4)
        ch(1,m2,k,1,1) = sn*(tr2+tr3)
        ch(1,m2,k,3,1) = sn*(tr2-tr3)
        ch(2,m2,k,1,1) = sn*(ti2+ti3)
        ch(2,m2,k,3,1) = sn*(ti2-ti3)
        ch(1,m2,k,2,1) = sn*(tr1+tr4)
        ch(1,m2,k,4,1) = sn*(tr1-tr4)
        ch(2,m2,k,2,1) = sn*(ti1+ti4)
        ch(2,m2,k,4,1) = sn*(ti1-ti4)
      end do
    end do

  else

    sn = 1.0E+00 / real ( 4 * l1 )

    do k = 1, l1
      do m1 = 1, m1d, im1
        ti1 = cc(2,m1,k,1,1)-cc(2,m1,k,1,3)
        ti2 = cc(2,m1,k,1,1)+cc(2,m1,k,1,3)
        tr4 = cc(2,m1,k,1,2)-cc(2,m1,k,1,4)
        ti3 = cc(2,m1,k,1,2)+cc(2,m1,k,1,4)
        tr1 = cc(1,m1,k,1,1)-cc(1,m1,k,1,3)
        tr2 = cc(1,m1,k,1,1)+cc(1,m1,k,1,3)
        ti4 = cc(1,m1,k,1,4)-cc(1,m1,k,1,2)
        tr3 = cc(1,m1,k,1,2)+cc(1,m1,k,1,4)
        cc(1,m1,k,1,1) = sn*(tr2+tr3)
        cc(1,m1,k,1,3) = sn*(tr2-tr3)
        cc(2,m1,k,1,1) = sn*(ti2+ti3)
        cc(2,m1,k,1,3) = sn*(ti2-ti3)
        cc(1,m1,k,1,2) = sn*(tr1+tr4)
        cc(1,m1,k,1,4) = sn*(tr1-tr4)
        cc(2,m1,k,1,2) = sn*(ti1+ti4)
        cc(2,m1,k,1,4) = sn*(ti1-ti4)
      end do
    end do

  end if

  return
end
subroutine cmf5kb ( lot, ido, l1, na, cc, im1, in1, ch, im2, in2, wa )

!*****************************************************************************80
!
!! CMF5KB is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(2,in1,l1,ido,5)
  real ch(2,in2,l1,5,ido)
  real chold1
  real chold2
  real ci2
  real ci3
  real ci4
  real ci5
  real cr2
  real cr3
  real cr4
  real cr5
  real di2
  real di3
  real di4
  real di5
  real dr2
  real dr3
  real dr4
  real dr5
  integer i
  integer im1
  integer im2
  integer k
  integer lot
  integer m1
  integer m1d
  integer m2
  integer m2s
  integer na
  real ti2
  real ti3
  real ti4
  real ti5
  real, parameter :: ti11 =  0.9510565162951536E+00
  real, parameter :: ti12 =  0.5877852522924731E+00
  real tr2
  real tr3
  real tr4
  real tr5
  real, parameter :: tr11 =  0.3090169943749474E+00
  real, parameter :: tr12 = -0.8090169943749474E+00
  real wa(ido,4,2)

  m1d = ( lot - 1 ) * im1 + 1
  m2s = 1 - im2

  if ( 1 < ido .or. na == 1 ) then

    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ti5 = cc(2,m1,k,1,2)-cc(2,m1,k,1,5)
        ti2 = cc(2,m1,k,1,2)+cc(2,m1,k,1,5)
        ti4 = cc(2,m1,k,1,3)-cc(2,m1,k,1,4)
        ti3 = cc(2,m1,k,1,3)+cc(2,m1,k,1,4)
        tr5 = cc(1,m1,k,1,2)-cc(1,m1,k,1,5)
        tr2 = cc(1,m1,k,1,2)+cc(1,m1,k,1,5)
        tr4 = cc(1,m1,k,1,3)-cc(1,m1,k,1,4)
        tr3 = cc(1,m1,k,1,3)+cc(1,m1,k,1,4)
        ch(1,m2,k,1,1) = cc(1,m1,k,1,1)+tr2+tr3
        ch(2,m2,k,1,1) = cc(2,m1,k,1,1)+ti2+ti3
        cr2 = cc(1,m1,k,1,1)+tr11*tr2+tr12*tr3
        ci2 = cc(2,m1,k,1,1)+tr11*ti2+tr12*ti3
        cr3 = cc(1,m1,k,1,1)+tr12*tr2+tr11*tr3
        ci3 = cc(2,m1,k,1,1)+tr12*ti2+tr11*ti3
        cr5 = ti11*tr5+ti12*tr4
        ci5 = ti11*ti5+ti12*ti4
        cr4 = ti12*tr5-ti11*tr4
        ci4 = ti12*ti5-ti11*ti4
        ch(1,m2,k,2,1) = cr2-ci5
        ch(1,m2,k,5,1) = cr2+ci5
        ch(2,m2,k,2,1) = ci2+cr5
        ch(2,m2,k,3,1) = ci3+cr4
        ch(1,m2,k,3,1) = cr3-ci4
        ch(1,m2,k,4,1) = cr3+ci4
        ch(2,m2,k,4,1) = ci3-cr4
        ch(2,m2,k,5,1) = ci2-cr5
      end do
    end do

    do i = 2, ido
      do k = 1, l1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ti5 = cc(2,m1,k,i,2)-cc(2,m1,k,i,5)
          ti2 = cc(2,m1,k,i,2)+cc(2,m1,k,i,5)
          ti4 = cc(2,m1,k,i,3)-cc(2,m1,k,i,4)
          ti3 = cc(2,m1,k,i,3)+cc(2,m1,k,i,4)
          tr5 = cc(1,m1,k,i,2)-cc(1,m1,k,i,5)
          tr2 = cc(1,m1,k,i,2)+cc(1,m1,k,i,5)
          tr4 = cc(1,m1,k,i,3)-cc(1,m1,k,i,4)
          tr3 = cc(1,m1,k,i,3)+cc(1,m1,k,i,4)
          ch(1,m2,k,1,i) = cc(1,m1,k,i,1)+tr2+tr3
          ch(2,m2,k,1,i) = cc(2,m1,k,i,1)+ti2+ti3
          cr2 = cc(1,m1,k,i,1)+tr11*tr2+tr12*tr3
          ci2 = cc(2,m1,k,i,1)+tr11*ti2+tr12*ti3
          cr3 = cc(1,m1,k,i,1)+tr12*tr2+tr11*tr3
          ci3 = cc(2,m1,k,i,1)+tr12*ti2+tr11*ti3
          cr5 = ti11*tr5+ti12*tr4
          ci5 = ti11*ti5+ti12*ti4
          cr4 = ti12*tr5-ti11*tr4
          ci4 = ti12*ti5-ti11*ti4
          dr3 = cr3-ci4
          dr4 = cr3+ci4
          di3 = ci3+cr4
          di4 = ci3-cr4
          dr5 = cr2+ci5
          dr2 = cr2-ci5
          di5 = ci2-cr5
          di2 = ci2+cr5
          ch(1,m2,k,2,i) = wa(i,1,1) * dr2 - wa(i,1,2) * di2
          ch(2,m2,k,2,i) = wa(i,1,1) * di2 + wa(i,1,2) * dr2
          ch(1,m2,k,3,i) = wa(i,2,1) * dr3 - wa(i,2,2) * di3
          ch(2,m2,k,3,i) = wa(i,2,1) * di3 + wa(i,2,2) * dr3
          ch(1,m2,k,4,i) = wa(i,3,1) * dr4 - wa(i,3,2) * di4
          ch(2,m2,k,4,i) = wa(i,3,1) * di4 + wa(i,3,2) * dr4
          ch(1,m2,k,5,i) = wa(i,4,1) * dr5 - wa(i,4,2) * di5
          ch(2,m2,k,5,i) = wa(i,4,1) * di5 + wa(i,4,2) * dr5
        end do
      end do
    end do

  else

    do k = 1, l1
      do m1 = 1, m1d, im1
        ti5 = cc(2,m1,k,1,2)-cc(2,m1,k,1,5)
        ti2 = cc(2,m1,k,1,2)+cc(2,m1,k,1,5)
        ti4 = cc(2,m1,k,1,3)-cc(2,m1,k,1,4)
        ti3 = cc(2,m1,k,1,3)+cc(2,m1,k,1,4)
        tr5 = cc(1,m1,k,1,2)-cc(1,m1,k,1,5)
        tr2 = cc(1,m1,k,1,2)+cc(1,m1,k,1,5)
        tr4 = cc(1,m1,k,1,3)-cc(1,m1,k,1,4)
        tr3 = cc(1,m1,k,1,3)+cc(1,m1,k,1,4)

        chold1 = cc(1,m1,k,1,1) + tr2 + tr3
        chold2 = cc(2,m1,k,1,1) + ti2 + ti3

        cr2 = cc(1,m1,k,1,1) + tr11 * tr2 + tr12 * tr3
        ci2 = cc(2,m1,k,1,1) + tr11 * ti2 + tr12 * ti3
        cr3 = cc(1,m1,k,1,1) + tr12 * tr2 + tr11 * tr3
        ci3 = cc(2,m1,k,1,1) + tr12 * ti2 + tr11 * ti3

        cc(1,m1,k,1,1) = chold1
        cc(2,m1,k,1,1) = chold2

        cr5 = ti11*tr5 + ti12*tr4
        ci5 = ti11*ti5 + ti12*ti4
        cr4 = ti12*tr5 - ti11*tr4
        ci4 = ti12*ti5 - ti11*ti4
        cc(1,m1,k,1,2) = cr2-ci5
        cc(1,m1,k,1,5) = cr2+ci5
        cc(2,m1,k,1,2) = ci2+cr5
        cc(2,m1,k,1,3) = ci3+cr4
        cc(1,m1,k,1,3) = cr3-ci4
        cc(1,m1,k,1,4) = cr3+ci4
        cc(2,m1,k,1,4) = ci3-cr4
        cc(2,m1,k,1,5) = ci2-cr5
      end do
    end do

  end if

  return
end
subroutine cmf5kf ( lot, ido, l1, na, cc, im1, in1, ch, im2, in2, wa )

!*****************************************************************************80
!
!! CMF5KF is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer l1

  real cc(2,in1,l1,ido,5)
  real ch(2,in2,l1,5,ido)
  real chold1
  real chold2
  real ci2
  real ci3
  real ci4
  real ci5
  real cr2
  real cr3
  real cr4
  real cr5
  real di2
  real di3
  real di4
  real di5
  real dr2
  real dr3
  real dr4
  real dr5
  integer i
  integer im1
  integer im2
  integer k
  integer lot
  integer m1
  integer m1d
  integer m2
  integer m2s
  integer na
  real sn
  real ti2
  real ti3
  real ti4
  real ti5
  real, parameter :: ti11 = -0.9510565162951536E+00
  real, parameter :: ti12 = -0.5877852522924731E+00
  real tr2
  real tr3
  real tr4
  real tr5
  real, parameter :: tr11 =  0.3090169943749474E+00
  real, parameter :: tr12 = -0.8090169943749474E+00
  real wa(ido,4,2)

  m1d = ( lot - 1 ) * im1 + 1
  m2s = 1 - im2

  if ( 1 < ido ) then

    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ti5 = cc(2,m1,k,1,2)-cc(2,m1,k,1,5)
        ti2 = cc(2,m1,k,1,2)+cc(2,m1,k,1,5)
        ti4 = cc(2,m1,k,1,3)-cc(2,m1,k,1,4)
        ti3 = cc(2,m1,k,1,3)+cc(2,m1,k,1,4)
        tr5 = cc(1,m1,k,1,2)-cc(1,m1,k,1,5)
        tr2 = cc(1,m1,k,1,2)+cc(1,m1,k,1,5)
        tr4 = cc(1,m1,k,1,3)-cc(1,m1,k,1,4)
        tr3 = cc(1,m1,k,1,3)+cc(1,m1,k,1,4)
        ch(1,m2,k,1,1) = cc(1,m1,k,1,1)+tr2+tr3
        ch(2,m2,k,1,1) = cc(2,m1,k,1,1)+ti2+ti3
        cr2 = cc(1,m1,k,1,1)+tr11*tr2+tr12*tr3
        ci2 = cc(2,m1,k,1,1)+tr11*ti2+tr12*ti3
        cr3 = cc(1,m1,k,1,1)+tr12*tr2+tr11*tr3
        ci3 = cc(2,m1,k,1,1)+tr12*ti2+tr11*ti3
        cr5 = ti11*tr5+ti12*tr4
        ci5 = ti11*ti5+ti12*ti4
        cr4 = ti12*tr5-ti11*tr4
        ci4 = ti12*ti5-ti11*ti4
        ch(1,m2,k,2,1) = cr2-ci5
        ch(1,m2,k,5,1) = cr2+ci5
        ch(2,m2,k,2,1) = ci2+cr5
        ch(2,m2,k,3,1) = ci3+cr4
        ch(1,m2,k,3,1) = cr3-ci4
        ch(1,m2,k,4,1) = cr3+ci4
        ch(2,m2,k,4,1) = ci3-cr4
        ch(2,m2,k,5,1) = ci2-cr5
      end do
    end do

    do i = 2, ido
      do k = 1, l1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ti5 = cc(2,m1,k,i,2)-cc(2,m1,k,i,5)
          ti2 = cc(2,m1,k,i,2)+cc(2,m1,k,i,5)
          ti4 = cc(2,m1,k,i,3)-cc(2,m1,k,i,4)
          ti3 = cc(2,m1,k,i,3)+cc(2,m1,k,i,4)
          tr5 = cc(1,m1,k,i,2)-cc(1,m1,k,i,5)
          tr2 = cc(1,m1,k,i,2)+cc(1,m1,k,i,5)
          tr4 = cc(1,m1,k,i,3)-cc(1,m1,k,i,4)
          tr3 = cc(1,m1,k,i,3)+cc(1,m1,k,i,4)
          ch(1,m2,k,1,i) = cc(1,m1,k,i,1)+tr2+tr3
          ch(2,m2,k,1,i) = cc(2,m1,k,i,1)+ti2+ti3
          cr2 = cc(1,m1,k,i,1)+tr11*tr2+tr12*tr3
          ci2 = cc(2,m1,k,i,1)+tr11*ti2+tr12*ti3
          cr3 = cc(1,m1,k,i,1)+tr12*tr2+tr11*tr3
          ci3 = cc(2,m1,k,i,1)+tr12*ti2+tr11*ti3
          cr5 = ti11*tr5+ti12*tr4
          ci5 = ti11*ti5+ti12*ti4
          cr4 = ti12*tr5-ti11*tr4
          ci4 = ti12*ti5-ti11*ti4
          dr3 = cr3-ci4
          dr4 = cr3+ci4
          di3 = ci3+cr4
          di4 = ci3-cr4
          dr5 = cr2+ci5
          dr2 = cr2-ci5
          di5 = ci2-cr5
          di2 = ci2+cr5
          ch(1,m2,k,2,i) = wa(i,1,1)*dr2+wa(i,1,2)*di2
          ch(2,m2,k,2,i) = wa(i,1,1)*di2-wa(i,1,2)*dr2
          ch(1,m2,k,3,i) = wa(i,2,1)*dr3+wa(i,2,2)*di3
          ch(2,m2,k,3,i) = wa(i,2,1)*di3-wa(i,2,2)*dr3
          ch(1,m2,k,4,i) = wa(i,3,1)*dr4+wa(i,3,2)*di4
          ch(2,m2,k,4,i) = wa(i,3,1)*di4-wa(i,3,2)*dr4
          ch(1,m2,k,5,i) = wa(i,4,1)*dr5+wa(i,4,2)*di5
          ch(2,m2,k,5,i) = wa(i,4,1)*di5-wa(i,4,2)*dr5
        end do
      end do
    end do

  else if ( na == 1 ) then

    sn = 1.0E+00 / real ( 5 * l1 )
 
    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ti5 = cc(2,m1,k,1,2)-cc(2,m1,k,1,5)
        ti2 = cc(2,m1,k,1,2)+cc(2,m1,k,1,5)
        ti4 = cc(2,m1,k,1,3)-cc(2,m1,k,1,4)
        ti3 = cc(2,m1,k,1,3)+cc(2,m1,k,1,4)
        tr5 = cc(1,m1,k,1,2)-cc(1,m1,k,1,5)
        tr2 = cc(1,m1,k,1,2)+cc(1,m1,k,1,5)
        tr4 = cc(1,m1,k,1,3)-cc(1,m1,k,1,4)
        tr3 = cc(1,m1,k,1,3)+cc(1,m1,k,1,4)
        ch(1,m2,k,1,1) = sn*(cc(1,m1,k,1,1)+tr2+tr3)
        ch(2,m2,k,1,1) = sn*(cc(2,m1,k,1,1)+ti2+ti3)
        cr2 = cc(1,m1,k,1,1)+tr11*tr2+tr12*tr3
        ci2 = cc(2,m1,k,1,1)+tr11*ti2+tr12*ti3
        cr3 = cc(1,m1,k,1,1)+tr12*tr2+tr11*tr3
        ci3 = cc(2,m1,k,1,1)+tr12*ti2+tr11*ti3
        cr5 = ti11*tr5+ti12*tr4
        ci5 = ti11*ti5+ti12*ti4
        cr4 = ti12*tr5-ti11*tr4
        ci4 = ti12*ti5-ti11*ti4
        ch(1,m2,k,2,1) = sn*(cr2-ci5)
        ch(1,m2,k,5,1) = sn*(cr2+ci5)
        ch(2,m2,k,2,1) = sn*(ci2+cr5)
        ch(2,m2,k,3,1) = sn*(ci3+cr4)
        ch(1,m2,k,3,1) = sn*(cr3-ci4)
        ch(1,m2,k,4,1) = sn*(cr3+ci4)
        ch(2,m2,k,4,1) = sn*(ci3-cr4)
        ch(2,m2,k,5,1) = sn*(ci2-cr5)
      end do
    end do

  else

    sn = 1.0E+00 / real ( 5 * l1 )

    do k = 1, l1
      do m1 = 1, m1d, im1

        ti5 = cc(2,m1,k,1,2) - cc(2,m1,k,1,5)
        ti2 = cc(2,m1,k,1,2) + cc(2,m1,k,1,5)
        ti4 = cc(2,m1,k,1,3) - cc(2,m1,k,1,4)
        ti3 = cc(2,m1,k,1,3) + cc(2,m1,k,1,4)
        tr5 = cc(1,m1,k,1,2) - cc(1,m1,k,1,5)
        tr2 = cc(1,m1,k,1,2) + cc(1,m1,k,1,5)
        tr4 = cc(1,m1,k,1,3) - cc(1,m1,k,1,4)
        tr3 = cc(1,m1,k,1,3) + cc(1,m1,k,1,4)

        chold1 = sn * ( cc(1,m1,k,1,1) + tr2 + tr3 )
        chold2 = sn * ( cc(2,m1,k,1,1) + ti2 + ti3 )

        cr2 = cc(1,m1,k,1,1) + tr11 * tr2 + tr12 * tr3
        ci2 = cc(2,m1,k,1,1) + tr11 * ti2 + tr12 * ti3
        cr3 = cc(1,m1,k,1,1) + tr12 * tr2 + tr11 * tr3
        ci3 = cc(2,m1,k,1,1) + tr12 * ti2 + tr11 * ti3

        cc(1,m1,k,1,1) = chold1
        cc(2,m1,k,1,1) = chold2

        cr5 = ti11 * tr5 + ti12 * tr4
        ci5 = ti11 * ti5 + ti12 * ti4
        cr4 = ti12 * tr5 - ti11 * tr4
        ci4 = ti12 * ti5 - ti11 * ti4

        cc(1,m1,k,1,2) = sn * ( cr2 - ci5 )
        cc(1,m1,k,1,5) = sn * ( cr2 + ci5 )
        cc(2,m1,k,1,2) = sn * ( ci2 + cr5 )
        cc(2,m1,k,1,3) = sn * ( ci3 + cr4 )
        cc(1,m1,k,1,3) = sn * ( cr3 - ci4 )
        cc(1,m1,k,1,4) = sn * ( cr3 + ci4 )
        cc(2,m1,k,1,4) = sn * ( ci3 - cr4 )
        cc(2,m1,k,1,5) = sn * ( ci2 - cr5 )

      end do
    end do

  end if

  return
end
subroutine cmfgkb ( lot, ido, ip, l1, lid, na, cc, cc1, im1, in1, &
  ch, ch1, im2, in2, wa )

!*****************************************************************************80
!
!! CMFGKB is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer ip
  integer l1
  integer lid

  real cc(2,in1,l1,ip,ido)
  real cc1(2,in1,lid,ip)
  real ch(2,in2,l1,ido,ip)
  real ch1(2,in2,lid,ip)
  real chold1
  real chold2
  integer i
  integer idlj
  integer im1
  integer im2
  integer ipp2
  integer ipph
  integer j
  integer jc
  integer k
  integer ki
  integer l
  integer lc
  integer lot
  integer m1
  integer m1d
  integer m2
  integer m2s
  integer na
  real wa(ido,ip-1,2)
  real wai
  real war

  m1d = ( lot - 1 ) * im1 + 1
  m2s = 1 - im2
  ipp2 = ip + 2
  ipph = ( ip + 1 ) / 2

  do ki = 1, lid
    m2 = m2s
    do m1 = 1, m1d, im1
      m2 = m2 + im2
      ch1(1,m2,ki,1) = cc1(1,m1,ki,1)
      ch1(2,m2,ki,1) = cc1(2,m1,ki,1)
    end do
  end do

  do j = 2, ipph
    jc = ipp2 - j
    do ki = 1, lid
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch1(1,m2,ki,j) =  cc1(1,m1,ki,j) + cc1(1,m1,ki,jc)
        ch1(1,m2,ki,jc) = cc1(1,m1,ki,j) - cc1(1,m1,ki,jc)
        ch1(2,m2,ki,j) =  cc1(2,m1,ki,j) + cc1(2,m1,ki,jc)
        ch1(2,m2,ki,jc) = cc1(2,m1,ki,j) - cc1(2,m1,ki,jc)
      end do
    end do
  end do

  do j = 2, ipph
    do ki = 1, lid
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        cc1(1,m1,ki,1) = cc1(1,m1,ki,1) + ch1(1,m2,ki,j)
        cc1(2,m1,ki,1) = cc1(2,m1,ki,1) + ch1(2,m2,ki,j)
      end do
    end do
  end do

  do l = 2, ipph

    lc = ipp2 - l
    do ki = 1, lid
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2

        cc1(1,m1,ki,l)  = ch1(1,m2,ki,1) + wa(1,l-1,1) * ch1(1,m2,ki,2)
        cc1(1,m1,ki,lc) =                  wa(1,l-1,2) * ch1(1,m2,ki,ip)
        cc1(2,m1,ki,l)  = ch1(2,m2,ki,1) + wa(1,l-1,1) * ch1(2,m2,ki,2)
        cc1(2,m1,ki,lc) =                  wa(1,l-1,2) * ch1(2,m2,ki,ip)

      end do
    end do

    do j = 3, ipph
      jc = ipp2 - j
      idlj = mod ( ( l - 1 ) * ( j - 1 ), ip )
      war = wa(1,idlj,1)
      wai = wa(1,idlj,2)
      do ki = 1, lid
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          cc1(1,m1,ki,l)  = cc1(1,m1,ki,l)  + war * ch1(1,m2,ki,j)
          cc1(1,m1,ki,lc) = cc1(1,m1,ki,lc) + wai * ch1(1,m2,ki,jc)
          cc1(2,m1,ki,l)  = cc1(2,m1,ki,l)  + war * ch1(2,m2,ki,j)
          cc1(2,m1,ki,lc) = cc1(2,m1,ki,lc) + wai * ch1(2,m2,ki,jc)
        end do
      end do
    end do

  end do

  if( 1 < ido .or. na == 1 ) then

    do ki = 1, lid
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch1(1,m2,ki,1) = cc1(1,m1,ki,1)
        ch1(2,m2,ki,1) = cc1(2,m1,ki,1)
      end do
    end do

    do j = 2, ipph
      jc = ipp2 - j
      do ki = 1, lid
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ch1(1,m2,ki,j)  = cc1(1,m1,ki,j) - cc1(2,m1,ki,jc)
          ch1(1,m2,ki,jc) = cc1(1,m1,ki,j) + cc1(2,m1,ki,jc)
          ch1(2,m2,ki,jc) = cc1(2,m1,ki,j) - cc1(1,m1,ki,jc)
          ch1(2,m2,ki,j)  = cc1(2,m1,ki,j) + cc1(1,m1,ki,jc)
        end do
      end do
    end do

    if ( ido == 1 ) then
      return
    end if

    do i = 1, ido
      do k = 1, l1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          cc(1,m1,k,1,i) = ch(1,m2,k,i,1)
          cc(2,m1,k,1,i) = ch(2,m2,k,i,1)
        end do
      end do
    end do

    do j = 2, ip
      do k = 1, l1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          cc(1,m1,k,j,1) = ch(1,m2,k,1,j)
          cc(2,m1,k,j,1) = ch(2,m2,k,1,j)
        end do
      end do
    end do

    do j = 2, ip
      do i = 2, ido
        do k = 1, l1
          m2 = m2s
          do m1 = 1, m1d, im1
            m2 = m2 + im2
            cc(1,m1,k,j,i) = wa(i,j-1,1) * ch(1,m2,k,i,j) &
                           - wa(i,j-1,2) * ch(2,m2,k,i,j)
            cc(2,m1,k,j,i) = wa(i,j-1,1) * ch(2,m2,k,i,j) &
                           + wa(i,j-1,2) * ch(1,m2,k,i,j)
          end do
        end do
      end do
    end do

  else

    do j = 2, ipph
      jc = ipp2 - j
      do ki = 1, lid
        do m1 = 1, m1d, im1

          chold1         = cc1(1,m1,ki,j) - cc1(2,m1,ki,jc)
          chold2         = cc1(1,m1,ki,j) + cc1(2,m1,ki,jc)
          cc1(1,m1,ki,j) = chold1

          cc1(2,m1,ki,jc) = cc1(2,m1,ki,j) - cc1(1,m1,ki,jc)
          cc1(2,m1,ki,j)  = cc1(2,m1,ki,j) + cc1(1,m1,ki,jc)
          cc1(1,m1,ki,jc) = chold2

        end do
      end do
    end do

  end if

  return
end
subroutine cmfgkf ( lot, ido, ip, l1, lid, na, cc, cc1, im1, in1, &
  ch, ch1, im2, in2, wa )

!*****************************************************************************80
!
!! CMFGKF is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer ido
  integer in1
  integer in2
  integer ip
  integer l1
  integer lid

  real cc(2,in1,l1,ip,ido)
  real cc1(2,in1,lid,ip)
  real ch(2,in2,l1,ido,ip)
  real ch1(2,in2,lid,ip)
  real chold1
  real chold2
  integer i
  integer idlj
  integer im1
  integer im2
  integer ipp2
  integer ipph
  integer j 
  integer jc
  integer k
  integer ki
  integer l
  integer lc
  integer lot
  integer m1
  integer m1d
  integer m2
  integer m2s
  integer na
  real sn
  real wa(ido,ip-1,2)
  real wai
  real war

  m1d = ( lot - 1 ) * im1 + 1
  m2s = 1 - im2
  ipp2 = ip + 2
  ipph = ( ip + 1 ) / 2

  do ki = 1, lid
    m2 = m2s
    do m1 = 1, m1d, im1
      m2 = m2 + im2
      ch1(1,m2,ki,1) = cc1(1,m1,ki,1)
      ch1(2,m2,ki,1) = cc1(2,m1,ki,1)
    end do
  end do

  do j = 2, ipph
    jc = ipp2 - j
    do ki = 1, lid
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch1(1,m2,ki,j) =  cc1(1,m1,ki,j) + cc1(1,m1,ki,jc)
        ch1(1,m2,ki,jc) = cc1(1,m1,ki,j) - cc1(1,m1,ki,jc)
        ch1(2,m2,ki,j) =  cc1(2,m1,ki,j) + cc1(2,m1,ki,jc)
        ch1(2,m2,ki,jc) = cc1(2,m1,ki,j) - cc1(2,m1,ki,jc)
      end do
    end do
  end do

  do j = 2, ipph
    do ki = 1, lid
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        cc1(1,m1,ki,1) = cc1(1,m1,ki,1) + ch1(1,m2,ki,j)
        cc1(2,m1,ki,1) = cc1(2,m1,ki,1) + ch1(2,m2,ki,j)
      end do
    end do
  end do

  do l = 2, ipph

    lc = ipp2 - l

    do ki = 1, lid
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        cc1(1,m1,ki,l)  = ch1(1,m2,ki,1) + wa(1,l-1,1) * ch1(1,m2,ki,2)
        cc1(1,m1,ki,lc) =                - wa(1,l-1,2) * ch1(1,m2,ki,ip)
        cc1(2,m1,ki,l)  = ch1(2,m2,ki,1) + wa(1,l-1,1) * ch1(2,m2,ki,2)
        cc1(2,m1,ki,lc) =                - wa(1,l-1,2) * ch1(2,m2,ki,ip)
      end do
    end do

    do j = 3, ipph
      jc = ipp2 - j
      idlj = mod ( ( l - 1 ) * ( j - 1 ), ip )
      war = wa(1,idlj,1)
      wai = -wa(1,idlj,2)
      do ki = 1, lid
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          cc1(1,m1,ki,l)  = cc1(1,m1,ki,l)  + war * ch1(1,m2,ki,j)
          cc1(1,m1,ki,lc) = cc1(1,m1,ki,lc) + wai * ch1(1,m2,ki,jc)
          cc1(2,m1,ki,l)  = cc1(2,m1,ki,l)  + war * ch1(2,m2,ki,j)
          cc1(2,m1,ki,lc) = cc1(2,m1,ki,lc) + wai * ch1(2,m2,ki,jc)
        end do
      end do
    end do

  end do

  if ( 1 < ido ) then

    do ki = 1, lid
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch1(1,m2,ki,1) = cc1(1,m1,ki,1)
        ch1(2,m2,ki,1) = cc1(2,m1,ki,1)
      end do
    end do

    do j = 2, ipph
      jc = ipp2 - j
      do ki = 1, lid
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ch1(1,m2,ki,j)  = cc1(1,m1,ki,j) - cc1(2,m1,ki,jc)
          ch1(2,m2,ki,j)  = cc1(2,m1,ki,j) + cc1(1,m1,ki,jc)
          ch1(1,m2,ki,jc) = cc1(1,m1,ki,j) + cc1(2,m1,ki,jc)
          ch1(2,m2,ki,jc) = cc1(2,m1,ki,j) - cc1(1,m1,ki,jc)
        end do
      end do
    end do

    do i = 1, ido
      do k = 1, l1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          cc(1,m1,k,1,i) = ch(1,m2,k,i,1)
          cc(2,m1,k,1,i) = ch(2,m2,k,i,1)
        end do
      end do
    end do

    do j = 2, ip
      do k = 1, l1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          cc(1,m1,k,j,1) = ch(1,m2,k,1,j)
          cc(2,m1,k,j,1) = ch(2,m2,k,1,j)
        end do
      end do
    end do

    do j = 2, ip
      do i = 2, ido
        do k = 1, l1
          m2 = m2s
          do m1 = 1, m1d, im1
            m2 = m2 + im2
            cc(1,m1,k,j,i) = wa(i,j-1,1) * ch(1,m2,k,i,j) &
                           + wa(i,j-1,2) * ch(2,m2,k,i,j)
            cc(2,m1,k,j,i) = wa(i,j-1,1) * ch(2,m2,k,i,j) &
                           - wa(i,j-1,2) * ch(1,m2,k,i,j)
          end do
        end do
      end do
    end do

  else if ( na == 1 ) then

    sn = 1.0E+00 / real ( ip * l1 )

    do ki = 1, lid
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch1(1,m2,ki,1) = sn * cc1(1,m1,ki,1)
        ch1(2,m2,ki,1) = sn * cc1(2,m1,ki,1)
      end do
    end do

    do j = 2, ipph
      jc = ipp2 - j
      do ki = 1, lid
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ch1(1,m2,ki,j)  = sn * ( cc1(1,m1,ki,j) - cc1(2,m1,ki,jc) )
          ch1(2,m2,ki,j)  = sn * ( cc1(2,m1,ki,j) + cc1(1,m1,ki,jc) )
          ch1(1,m2,ki,jc) = sn * ( cc1(1,m1,ki,j) + cc1(2,m1,ki,jc) )
          ch1(2,m2,ki,jc) = sn * ( cc1(2,m1,ki,j) - cc1(1,m1,ki,jc) )
        end do
      end do
    end do

  else

    sn = 1.0E+00 / real ( ip * l1 )

    do ki = 1, lid
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        cc1(1,m1,ki,1) = sn * cc1(1,m1,ki,1)
        cc1(2,m1,ki,1) = sn * cc1(2,m1,ki,1)
      end do
    end do

    do j = 2, ipph
      jc = ipp2 - j
      do ki = 1, lid
        do m1 = 1, m1d, im1
          chold1 = sn * ( cc1(1,m1,ki,j) - cc1(2,m1,ki,jc) )
          chold2 = sn * ( cc1(1,m1,ki,j) + cc1(2,m1,ki,jc) )
          cc1(1,m1,ki,j) = chold1
          cc1(2,m1,ki,jc) = sn * ( cc1(2,m1,ki,j) - cc1(1,m1,ki,jc) )
          cc1(2,m1,ki,j)  = sn * ( cc1(2,m1,ki,j) + cc1(1,m1,ki,jc) )
          cc1(1,m1,ki,jc) = chold2
        end do
      end do
    end do

  end if

  return
end
subroutine cmfm1b ( lot, jump, n, inc, c, ch, wa, fnf, fac )

!*****************************************************************************80
!
!! CMFM1B is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  complex c(*)
  real ch(*)
  real fac(*)
  real fnf
  integer ido
  integer inc
  integer ip
  integer iw
  integer jump
  integer k1
  integer l1
  integer l2
  integer lid
  integer lot
  integer n
  integer na
  integer nbr
  integer nf
  real wa(*)

  nf = int ( fnf )
  na = 0
  l1 = 1
  iw = 1

  do k1 = 1, nf

    ip = int ( fac(k1) )
    l2 = ip * l1
    ido = n / l2
    lid = l1 * ido
    nbr = 1 + na + 2 * min ( ip - 2, 4 )

    if ( nbr == 1 ) then
      call cmf2kb ( lot, ido, l1, na, c, jump, inc, ch, 1, lot, wa(iw) )
    else if ( nbr == 2 ) then
      call cmf2kb ( lot, ido, l1, na, ch, 1, lot, c, jump, inc, wa(iw) )
    else if ( nbr == 3 ) then
      call cmf3kb ( lot, ido, l1, na, c, jump, inc, ch, 1, lot, wa(iw) )
    else if ( nbr == 4 ) then
      call cmf3kb ( lot, ido, l1, na, ch, 1, lot, c, jump, inc, wa(iw) )
    else if ( nbr == 5 ) then
      call cmf4kb ( lot, ido, l1, na, c, jump, inc, ch, 1, lot, wa(iw) )
    else if ( nbr == 6 ) then
      call cmf4kb ( lot, ido, l1, na, ch, 1, lot, c, jump, inc, wa(iw) )
    else if ( nbr == 7 ) then
      call cmf5kb ( lot, ido, l1, na, c, jump, inc, ch, 1, lot, wa(iw) )
    else if ( nbr == 8 ) then
      call cmf5kb ( lot, ido, l1, na, ch, 1, lot, c, jump, inc, wa(iw) )
    else if ( nbr == 9 ) then
      call cmfgkb ( lot, ido, ip, l1, lid, na, c, c, jump, inc, ch, ch, &
        1, lot, wa(iw) )
    else if ( nbr == 10 ) then
      call cmfgkb ( lot, ido, ip, l1, lid, na, ch, ch, 1, lot, c, c, &
        jump, inc, wa(iw) )
    end if

    l1 = l2
    iw = iw + ( ip - 1 ) * ( ido + ido )

    if ( ip <= 5 ) then
      na = 1 - na
    end if

  end do

  return
end
subroutine cmfm1f ( lot, jump, n, inc, c, ch, wa, fnf, fac )

!*****************************************************************************80
!
!! CMFM1F is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  complex c(*)
  real ch(*)
  real fac(*)
  real fnf
  integer ido
  integer inc
  integer ip
  integer iw
  integer jump
  integer k1
  integer l1
  integer l2
  integer lid
  integer lot
  integer n
  integer na
  integer nbr
  integer nf
  real wa(*)

  nf = int ( fnf )
  na = 0
  l1 = 1
  iw = 1

  do k1 = 1, nf

    ip = int ( fac(k1) )
    l2 = ip * l1
    ido = n / l2
    lid = l1 * ido
    nbr = 1 + na + 2 * min ( ip - 2, 4 )

    if ( nbr == 1 ) then
      call cmf2kf ( lot, ido, l1, na, c, jump, inc, ch, 1, lot, wa(iw) )
    else if ( nbr == 2 ) then
      call cmf2kf ( lot, ido, l1, na, ch, 1, lot, c, jump, inc, wa(iw) )
    else if ( nbr == 3 ) then
      call cmf3kf ( lot, ido, l1, na, c, jump, inc, ch, 1, lot, wa(iw) )
    else if ( nbr == 4 ) then
      call cmf3kf ( lot, ido, l1, na, ch, 1, lot, c, jump, inc, wa(iw) )
    else if ( nbr == 5 ) then
      call cmf4kf ( lot, ido, l1, na, c, jump, inc, ch, 1, lot, wa(iw) )
    else if ( nbr == 6 ) then
      call cmf4kf ( lot, ido, l1, na, ch, 1, lot, c, jump, inc, wa(iw) )
    else if ( nbr == 7 ) then
      call cmf5kf ( lot, ido, l1, na, c, jump, inc, ch, 1, lot, wa(iw) )
    else if ( nbr == 8 ) then
      call cmf5kf ( lot, ido, l1, na, ch, 1, lot, c, jump, inc, wa(iw) )
    else if ( nbr == 9 ) then
      call cmfgkf ( lot, ido, ip, l1, lid, na, c, c, jump, inc, ch, ch, &
        1, lot, wa(iw) )
    else if ( nbr == 10 ) then
      call cmfgkf ( lot, ido, ip, l1, lid, na, ch, ch, 1, lot, c, c, &
        jump, inc, wa(iw) )
    end if

    l1 = l2
    iw = iw + ( ip - 1 ) * ( ido + ido )

    if ( ip <= 5 ) then
      na = 1 - na
    end if

  end do

  return
end
subroutine cosq1b ( n, inc, x, lenx, wsave, lensav, work, lenwrk, ier )

!*****************************************************************************80
!
!! COSQ1B: real backward cosine quarter wave transform, 1D.
!
!  Discussion:
!
!    COSQ1B computes the one-dimensional Fourier transform of a sequence 
!    which is a cosine series with odd wave numbers.  This transform is 
!    referred to as the backward transform or Fourier synthesis, transforming
!    the sequence from spectral to physical space.
!
!    This transform is normalized since a call to COSQ1B followed
!    by a call to COSQ1F (or vice-versa) reproduces the original
!    array  within roundoff error.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Modified:
!
!    31 March 2005
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
!    Input, integer N, the number of elements to be transformed in the 
!    sequence.  The transform is most efficient when N is a 
!    product of small primes.
!
!    Input, integer INC, the increment between the locations, in array R,
!    of two consecutive elements within the sequence.
!
!    Input/output, real R(LENR); on input, containing the sequence to be 
!    transformed, and on output, containing the transformed sequence.
!
!    Input, integer LENR, the dimension of the R array.  LENR must be at least
!    INC*(N-1)+ 1.
!
!    Input, real WSAVE(LENSAV).  WSAVE's contents must be initialized 
!    with a call to COSQ1I before the first call to routine COSQ1F or 
!    COSQ1B for a given transform length N.  WSAVE's contents may be 
!    re-used for subsequent calls to COSQ1F and COSQ1B with the same N.
!
!    Input, integer LENSAV, the dimension of the WSAVE array.  LENSAV must 
!    be at least 2*N + INT(LOG(REAL(N))) + 4.
!
!    Workspace, real WORK(LENWRK).
!
!    Input, integer LENWRK, the dimension of the WORK array.  LENWRK must be at 
!    least N.
!
!    Output, integer IER, error flag.
!    0, successful exit;
!    1, input parameter LENR   not big enough;
!    2, input parameter LENSAV not big enough;
!    3, input parameter LENWRK not big enough;
!    20, input error returned by lower level routine.
!
  implicit none

  integer inc
  integer lensav
  integer lenwrk

  integer ier
  integer ier1
  integer lenx
  integer n
  real ssqrt2
  real work(lenwrk)
  real wsave(lensav)
  real x(inc,*)
  real x1

  ier = 0

  if ( lenx < inc * ( n - 1 ) + 1 ) then
    ier = 1
    call xerfft ( 'COSQ1B', 6 )
    return
  end if

  if ( lensav < 2 * n + int ( log ( real ( n ) ) ) + 4 ) then
    ier = 2
    call xerfft ( 'COSQ1B', 8 )
    return
  end if

  if ( lenwrk < n ) then
    ier = 3
    call xerfft ( 'COSQ1B', 10 ) 
    return
  end if

  if ( n < 2 ) then
    return
  end if

  if ( n == 2 ) then
    ssqrt2 = 1.0E+00 / sqrt ( 2.0E+00 )
    x1 = x(1,1) + x(1,2)
    x(1,2) = ssqrt2 * ( x(1,1) - x(1,2) )
    x(1,1) = x1
    return
  end if

  call cosqb1 ( n, inc, x, wsave, work, ier1 )

  if ( ier1 /= 0 ) then
    ier = 20
    call xerfft ( 'COSQ1B', -5 )
    return
  end if

  return
end
subroutine cosq1f ( n, inc, x, lenx, wsave, lensav, work, lenwrk, ier )

!*****************************************************************************80
!
!! COSQ1F: real forward cosine quarter wave transform, 1D.
!
!  Discussion:
!
!    COSQ1F computes the one-dimensional Fourier transform of a sequence 
!    which is a cosine series with odd wave numbers.  This transform is 
!    referred to as the forward transform or Fourier analysis, transforming 
!    the sequence from physical to spectral space.
!
!    This transform is normalized since a call to COSQ1F followed
!    by a call to COSQ1B (or vice-versa) reproduces the original
!    array  within roundoff error.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Modified:
!
!    31 March 2005
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
!    Input, integer N, the number of elements to be transformed in the 
!    sequence.  The transform is most efficient when N is a 
!    product of small primes.
!
!    Input, integer INC, the increment between the locations, in array R,
!    of two consecutive elements within the sequence.
!
!    Input/output, real R(LENR); on input, containing the sequence to be 
!    transformed, and on output, containing the transformed sequence.
!
!    Input, integer LENR, the dimension of the R array.  LENR must be at least
!    INC*(N-1)+ 1.
!
!    Input, real WSAVE(LENSAV).  WSAVE's contents must be initialized 
!    with a call to COSQ1I before the first call to routine COSQ1F or 
!    COSQ1B for a given transform length N.  WSAVE's contents may be 
!    re-used for subsequent calls to COSQ1F and COSQ1B with the same N.
!
!    Input, integer LENSAV, the dimension of the WSAVE array.  LENSAV must 
!    be at least 2*N + INT(LOG(REAL(N))) + 4.
!
!    Workspace, real WORK(LENWRK).
!
!    Input, integer LENWRK, the dimension of the WORK array.  LENWRK must be at 
!    least N.
!
!    Output, integer IER, error flag.
!    0, successful exit;
!    1, input parameter LENR   not big enough;
!    2, input parameter LENSAV not big enough;
!    3, input parameter LENWRK not big enough;
!    20, input error returned by lower level routine.
!
  implicit none

  integer inc
  integer lensav
  integer lenwrk

  integer ier
  integer ier1
  integer n
  integer lenx
  real ssqrt2
  real tsqx
  real work(lenwrk)
  real wsave(lensav)
  real x(inc,*)

  ier = 0

  if ( lenx < inc * ( n - 1 ) + 1 ) then
    ier = 1
    call xerfft ( 'cosq1f', 6 )
    return
  end if

  if ( lensav < 2 * n + int ( log ( real ( n ) ) ) + 4 ) then
    ier = 2
    call xerfft ( 'cosq1f', 8 )
    return
  end if

  if ( lenwrk < n ) then
    ier = 3
    call xerfft ( 'cosq1f', 10 )
    return
  end if

  if ( n < 2 ) then
    return
  end if

  if ( n == 2 ) then
    ssqrt2 = 1.0E+00 / sqrt ( 2.0E+00 )
    tsqx = ssqrt2 * x(1,2)
    x(1,2) = 0.5E+00 * x(1,1) - tsqx
    x(1,1) = 0.5E+00 * x(1,1) + tsqx
    return
  end if

  call cosqf1 ( n, inc, x, wsave, work, ier1 )

  if ( ier1 /= 0 ) then
    ier = 20
    call xerfft ( 'cosq1f', -5 )
    return
  end if

  return
end
subroutine cosq1i ( n, wsave, lensav, ier )

!*****************************************************************************80
!
!! COSQ1I: initialization for COSQ1B and COSQ1F.
!
!  Discussion:
!
!    COSQ1I initializes array WSAVE for use in its companion routines 
!    COSQ1F and COSQ1B.  The prime factorization of N together with a 
!    tabulation of the trigonometric functions are computed and stored 
!    in array WSAVE.  Separate WSAVE arrays are required for different 
!    values of N.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Modified:
!
!    31 March 2005
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
!    Input, integer N, the length of the sequence to be transformed.  The 
!    transform is most efficient when N is a product of small primes.
!
!    Input, integer LENSAV, the dimension of the WSAVE array.  LENSAV must 
!    be at least 2*N + INT(LOG(REAL(N))) + 4.
!
!    Output, real WSAVE(LENSAV), containing the prime factors of N and
!    also containing certain trigonometric values which will be used 
!    in routines COSQ1B or COSQ1F.
!
!    Output, integer IER, error flag.
!    0, successful exit;
!    2, input parameter LENSAV not big enough;
!    20, input error returned by lower level routine.
!
  implicit none

  integer lensav

  real dt
  real fk
  integer ier
  integer ier1
  integer k
  integer lnsv
  integer n
  real pih
  real wsave(lensav)

  ier = 0

  if ( lensav < 2 * n + int ( log ( real ( n ) ) ) + 4 ) then
    ier = 2
    call xerfft ( 'cosq1i', 3 )
    return
  end if

  pih = 2.0E+00 * atan ( 1.0E+00 )
  dt = pih / real ( n )
  fk = 0.0E+00

  do k = 1, n
    fk = fk + 1.0E+00
    wsave(k) = cos ( fk * dt )
  end do

  lnsv = n + int ( log ( real ( n ) ) ) + 4

  call rfft1i ( n, wsave(n+1), lnsv, ier1 )

  if ( ier1 /= 0 ) then
    ier = 20
    call xerfft ( 'cosq1i', -5 )
    return
  end if

  return
end
subroutine cosqb1 ( n, inc, x, wsave, work, ier )

!*****************************************************************************80
!
!! COSQB1 is an FFTPACK5 auxiliary routine.
!
!  License:
!
!    Licensed under the GNU General Public License (GPL).
!    Copyright (C) 1995-2004, Scientific Computing Division,
!    University Corporation for Atmospheric Research
!
!  Author:
!
!    Paul Swarztrauber
!    Richard Valent
!
!  Reference:
!
!    Paul Swarztrauber,
!    Vectorizing the Fast Fourier Transforms,
!    in Parallel Computations,
!    edited by G. Rodrigue,
!    Academic Press, 1982.
!
!    Paul Swarztrauber,
!    Fast Fourier Transform Algorithms for Vector Computers,
!    Parallel Computing, pages 45-63, 1984.
!
!  Parameters:
!
  implicit none

  integer inc

  integer i
  integer ier
  integer ier1
  integer k
  integer kc
  integer lenx
  integer lnsv
  integer lnwk
  integer modn
  integer n
  integer np2
  integer ns2
  real work(*)
  real wsave(*)
  real x(inc,*)
  real xim1

  ier = 0
  ns2 = ( n + 1 ) / 2
  np2 = n + 2

  do i = 3, n, 2
    xim1 = x(1,i-1) + x(1,i)
    x(1,i) = 0.5E+00 * ( x(1,i-1) - x(1,i) )
    x(1,i-1) = 0.5E+00 * xim1
  end do

  x(1,1) = 0.5E+00 * x(1,1)
  modn = mod ( n, 2 )

  if ( modn == 0 ) then
    x(1,n) = 0.5E+00 * x(1,n)
  end if

  lenx = inc * ( n - 1 )  + 1
  lnsv = n + int ( log ( real ( n ) ) ) + 4
  lnwk = n

  call rfft1b ( n, inc, x, lenx, wsave(n+1), lnsv, work, lnwk, ier1 )

  if ( ier1 /= 0 ) then
    ier = 20
    call xerfft ( 'cosqb1', -5 )
    return
  end if

  do k = 2, ns2
    kc = np2 - k
    work(k)  = wsave(k-1) * x(1,kc) + wsave(kc-1) * x(1,k)
    work(kc) = wsave(k-1) * x(1,k)  - wsave