module fftpack5
!-----------------------------------------------------------------------
! Purpose: 
! 
! selected Fast Fourier Transform routines taken from FFTPACK5.1
! Used for Explicit Scalar Momentum Tendency (ESMT) code 
!
! Revision history: 
! Jan, 2018 - Walter Hannah
!              initial version - adapted from https://michalakes.svn.cloudforge.com/rrtmmic/HWRF/external/fftpack/fftpack5
!
!---------------------------------------------------------------------------
  
  use params, only: crm_rknd

  implicit none

  public rfftmi
  public rfftmf
  public rfftmb

  contains
  
!========================================================================================
! RFFTMI: initialization for RFFTMB and RFFTMF.
!========================================================================================
subroutine rfftmi ( n, wsave, lensav, ier )

  !*****************************************************************************80
  !
  !  Discussion:
  !
  !    RFFTMI initializes array WSAVE for use in its companion routines 
  !    RFFTMB and RFFTMF.  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 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) N, the length of each sequence to be 
  !    transformed.  The transform is most efficient when N is a product of 
  !    small primes.
  !
  !    Input, integer ( kind = 4 ) LENSAV, the dimension of the WSAVE array.  
  !    LENSAV must be at least N + INT(LOG(REAL(N))) + 4.
  !
  !    Output, real ( kind = crm_rknd ) WSAVE(LENSAV), work array containing the prime 
  !    factors of N and also containing certain trigonometric 
  !    values which will be used in routines RFFTMB or RFFTMF.
  !
  !    Output, integer ( kind = 4 ) IER, error flag.
  !    0, successful exit;
  !    2, input parameter LENSAV not big enough.
  !
  implicit none

  integer ( kind = 4 )     ,intent(in   ) :: lensav
  integer ( kind = 4 )     ,intent(in   ) :: n
  integer ( kind = 4 )     ,intent(  out) :: ier
  real ( kind = crm_rknd ) ,intent(inout) :: wsave(lensav)

  ier = 0

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

  if ( n == 1 ) then
    return
  end if

  call mrfti1 ( n, wsave(1), wsave(n+1) )

  return
end
!========================================================================================
! RFFTMF: real single precision forward FFT, 1D, multiple vectors.
!========================================================================================
subroutine rfftmf ( lot, jump, n, inc, r, lenr, wsave, lensav, &
              work, lenwrk, ier )

  !*****************************************************************************80
  !
  !  Discussion:
  !
  !    RFFTMF computes the one-dimensional Fourier transform of multiple 
  !    periodic sequences within a real 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 RFFTMF followed
  !    by a call to RFFTMB (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 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) LOT, the number of sequences to be transformed 
  !    within array R.
  !
  !    Input, integer ( kind = 4 ) JUMP, the increment between the locations, in 
  !    array R, of the first elements of two consecutive sequences to be 
  !    transformed.
  !
  !    Input, integer ( kind = 4 ) N, the length of each sequence to be 
  !    transformed.  The transform is most efficient when N is a product of 
  !    small primes.
  !
  !    Input, integer ( kind = 4 ) INC, the increment between the locations, 
  !    in array R, of two consecutive elements within the same sequence.
  !
  !    Input/output, real ( kind = crm_rknd ) R(LENR), real array containing LOT 
  !    sequences, each having length N.  R can have any number of dimensions, but 
  !    the total number of locations must be at least LENR.  On input, the
  !    physical data to be transformed, on output the spectral data.
  !
  !    Input, integer ( kind = 4 ) LENR, the dimension of the R array.  
  !    LENR must be at least (LOT-1)*JUMP + INC*(N-1) + 1.
  !
  !    Input, real ( kind = crm_rknd ) WSAVE(LENSAV).  WSAVE's contents must be 
  !    initialized with a call to RFFTMI before the first call to routine RFFTMF 
  !    or RFFTMB for a given transform length N.  
  !
  !    Input, integer ( kind = 4 ) LENSAV, the dimension of the WSAVE array.  
  !    LENSAV must be at least N + INT(LOG(REAL(N))) + 4.
  !
  !    Workspace, real ( kind = crm_rknd ) WORK(LENWRK).
  !
  !    Input, integer ( kind = 4 ) LENWRK, the dimension of the WORK array.  
  !    LENWRK must be at least LOT*N.
  !
  !    Output, integer ( kind = 4 ) 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;
  !    4, input parameters INC, JUMP, N, LOT are not consistent.
  !
  implicit none
  
  integer ( kind = 4 )        ,intent(in   ) :: lot
  integer ( kind = 4 )        ,intent(in   ) :: jump
  integer ( kind = 4 )        ,intent(in   ) :: n
  integer ( kind = 4 )        ,intent(in   ) :: inc
  integer ( kind = 4 )        ,intent(in   ) :: lenr
  integer ( kind = 4 )        ,intent(in   ) :: lensav
  integer ( kind = 4 )        ,intent(in   ) :: lenwrk
  integer ( kind = 4 )        ,intent(  out) :: ier
  real    ( kind = crm_rknd ) ,intent(inout) :: r(lenr)
  real    ( kind = crm_rknd ) ,intent(inout) :: work(lenwrk)
  real    ( kind = crm_rknd ) ,intent(inout) :: wsave(lensav)

  ier = 0

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

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

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

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

  if ( n == 1 ) then
    return
  end if

  call mrftf1 ( lot, jump, n, inc, r, work, wsave, wsave(n+1) )

  return
end
!========================================================================================
! RFFTMB: real single precision backward FFT, 1D, multiple vectors.
!========================================================================================
subroutine rfftmb ( lot, jump, n, inc, r, lenr, wsave, lensav, &
              work, lenwrk, ier )

  !*****************************************************************************80
  !
  !  Discussion:
  !
  !    RFFTMB computes the one-dimensional Fourier transform of multiple 
  !    periodic sequences within a real 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 RFFTMB followed
  !    by a call to RFFTMF (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 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) LOT, the number of sequences to be transformed 
  !    within array R.
  !
  !    Input, integer ( kind = 4 ) JUMP, the increment between the locations, in
  !    array R, of the first elements of two consecutive sequences to be 
  !    transformed.
  !
  !    Input, integer ( kind = 4 ) N, the length of each sequence to be 
  !    transformed.  The transform is most efficient when N is a product of 
  !    small primes.
  !
  !    Input, integer ( kind = 4 ) INC, the increment between the locations, in 
  !    array R, of two consecutive elements within the same sequence.
  !
  !    Input/output, real ( kind = crm_rknd ) R(LENR), real array containing LOT 
  !    sequences, each having length N.  R can have any number of dimensions, 
  !    but the total number of locations must be at least LENR.  On input, the
  !    spectral data to be transformed, on output the physical data.
  !
  !    Input, integer ( kind = 4 ) LENR, the dimension of the R array. 
  !    LENR must be at least (LOT-1)*JUMP + INC*(N-1) + 1.
  !
  !    Input, real ( kind = crm_rknd ) WSAVE(LENSAV).  WSAVE's contents must be 
  !    initialized with a call to RFFTMI before the first call to routine RFFTMF 
  !    or RFFTMB for a given transform length N.  
  !
  !    Input, integer ( kind = 4 ) LENSAV, the dimension of the WSAVE array. 
  !    LENSAV must  be at least N + INT(LOG(REAL(N))) + 4.
  !
  !    Workspace, real ( kind = crm_rknd ) WORK(LENWRK).
  !
  !    Input, integer ( kind = 4 ) LENWRK, the dimension of the WORK array.  
  !    LENWRK must be at least LOT*N.
  !
  !    Output, integer ( kind = 4 ) 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;
  !    4, input parameters INC, JUMP, N, LOT are not consistent.
  !
  implicit none

  integer ( kind = 4 )        ,intent(in   ) :: lenr
  integer ( kind = 4 )        ,intent(in   ) :: lensav
  integer ( kind = 4 )        ,intent(in   ) :: lenwrk
  integer ( kind = 4 )        ,intent(in   ) :: inc
  integer ( kind = 4 )        ,intent(in   ) :: jump
  integer ( kind = 4 )        ,intent(in   ) :: lot
  integer ( kind = 4 )        ,intent(in   ) :: n
  integer ( kind = 4 )        ,intent(  out) :: ier
  real    ( kind = crm_rknd ) ,intent(inout) :: r(lenr)
  real    ( kind = crm_rknd ) ,intent(inout) :: work(lenwrk)
  real    ( kind = crm_rknd ) ,intent(inout) :: wsave(lensav)

  ier = 0

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

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

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

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

  if ( n == 1 ) then
    return
  end if

  call mrftb1 ( lot, jump, n, inc, r, work, wsave, wsave(n+1) )

  return
end

!========================================================================================
! MRADB2 is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mradb2 (m,ido,l1,cc,im1,in1,ch,im2,in2,wa1)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 )     ,intent(in   ) :: ido
  integer ( kind = 4 )     ,intent(in   ) :: in1
  integer ( kind = 4 )     ,intent(in   ) :: in2
  integer ( kind = 4 )     ,intent(in   ) :: l1
  integer ( kind = 4 )     ,intent(in   ) :: im1
  integer ( kind = 4 )     ,intent(in   ) :: im2
  integer ( kind = 4 )     ,intent(in   ) :: m
  real ( kind = crm_rknd ) ,intent(in   ) :: cc(in1,ido,2,l1)
  real ( kind = crm_rknd ) ,intent(in   ) :: wa1(ido)
  real ( kind = crm_rknd ) ,intent(inout) :: ch(in2,ido,l1,2)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) idp2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m1
  integer ( kind = 4 ) m1d
  integer ( kind = 4 ) m2
  integer ( kind = 4 ) m2s
  

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

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

  if ( ido < 2 ) then
    return
  end if

  if ( 2 < ido ) then

    idp2 = ido + 2

    do k = 1, l1
      do i = 3, ido, 2
        ic = idp2 - i
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ch(m2,i-1,k,1) = cc(m1,i-1,1,k) + cc(m1,ic-1,2,k)
          ch(m2,i,k,1)   = cc(m1,i,1,k) - cc(m1,ic,2,k)
          ch(m2,i-1,k,2) = wa1(i-2) * ( cc(m1,i-1,1,k) - cc(m1,ic-1,2,k) ) &
                         - wa1(i-1) * ( cc(m1,i,1,k)   + cc(m1,ic,2,k) )
          ch(m2,i,k,2)   = wa1(i-2) * ( cc(m1,i,1,k)   + cc(m1,ic,2,k) ) &
                         + wa1(i-1) * ( cc(m1,i-1,1,k) - cc(m1,ic-1,2,k) )
        end do
      end do
    end do

    if ( mod ( ido, 2 ) == 1 ) then
      return
    end if

  end if

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

  return
end
!========================================================================================
! MRADB3 is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mradb3 (m,ido,l1,cc,im1,in1,ch,im2,in2,wa1,wa2)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) ido
  integer ( kind = 4 ) in1
  integer ( kind = 4 ) in2
  integer ( kind = 4 ) l1

  real ( kind = crm_rknd ) arg
  real ( kind = crm_rknd ) ,intent(in   ) :: cc(in1,ido,3,l1)
  real ( kind = crm_rknd ) ,intent(inout) :: ch(in2,ido,l1,3)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) idp2
  integer ( kind = 4 ) im1
  integer ( kind = 4 ) im2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) m1
  integer ( kind = 4 ) m1d
  integer ( kind = 4 ) m2
  integer ( kind = 4 ) m2s
  real ( kind = crm_rknd ) taui
  real ( kind = crm_rknd ) taur
  real ( kind = crm_rknd ) wa1(ido)
  real ( kind = crm_rknd ) wa2(ido)

  m1d = ( m - 1 ) * im1 + 1
  m2s = 1 - im2
  arg = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 ) / 3.0E+00
  taur = cos ( arg )
  taui = sin ( arg )

  do k = 1, l1
    m2 = m2s
    do m1 = 1, m1d, im1
      m2 = m2 + im2
      ch(m2,1,k,1) = cc(m1,1,1,k) + 2.0E+00 * cc(m1,ido,2,k)
      ch(m2,1,k,2) = cc(m1,1,1,k) + ( 2.0E+00 * taur ) * cc(m1,ido,2,k) &
        - ( 2.0E+00 * taui ) * cc(m1,1,3,k)
      ch(m2,1,k,3) = cc(m1,1,1,k) + ( 2.0E+00 * taur ) * cc(m1,ido,2,k) &
        + 2.0E+00 * taui * cc(m1,1,3,k)
    end do
  end do

  if ( ido == 1 ) then
    return
  end if

  idp2 = ido + 2

  do k = 1, l1
    do i = 3, ido, 2
      ic = idp2 - i
      m2 = m2s
      do m1 = 1, m1d, im1

        m2 = m2 + im2

        ch(m2,i-1,k,1) = cc(m1,i-1,1,k)+(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k))

        ch(m2,i,k,1) = cc(m1,i,1,k)+(cc(m1,i,3,k)-cc(m1,ic,2,k))

        ch(m2,i-1,k,2) = wa1(i-2)* &
          ((cc(m1,i-1,1,k)+taur*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)))- &
          (taui*(cc(m1,i,3,k)+cc(m1,ic,2,k)))) - wa1(i-1)* &
          ((cc(m1,i,1,k)+taur*(cc(m1,i,3,k)-cc(m1,ic,2,k)))+ &
          (taui*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))))

        ch(m2,i,k,2) = wa1(i-2)* &
          ((cc(m1,i,1,k)+taur*(cc(m1,i,3,k)-cc(m1,ic,2,k)))+ &
          (taui*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k)))) + wa1(i-1)* &
          ((cc(m1,i-1,1,k)+taur*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)))- &
          (taui*(cc(m1,i,3,k)+cc(m1,ic,2,k))))

        ch(m2,i-1,k,3) = wa2(i-2)* &
          ((cc(m1,i-1,1,k)+taur*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)))+ &
          (taui*(cc(m1,i,3,k)+cc(m1,ic,2,k)))) - wa2(i-1)* &
          ((cc(m1,i,1,k)+taur*(cc(m1,i,3,k)-cc(m1,ic,2,k)))- &
          (taui*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))))

        ch(m2,i,k,3) = wa2(i-2)* &
          ((cc(m1,i,1,k)+taur*(cc(m1,i,3,k)-cc(m1,ic,2,k)))- &
          (taui*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k)))) + wa2(i-1)* &
          ((cc(m1,i-1,1,k)+taur*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)))+ &
          (taui*(cc(m1,i,3,k)+cc(m1,ic,2,k))))

      end do
    end do
  end do

  return
end
!========================================================================================
! MRADB4 is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mradb4 (m,ido,l1,cc,im1,in1,ch,im2,in2,wa1,wa2,wa3)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) ido
  integer ( kind = 4 ) in1
  integer ( kind = 4 ) in2
  integer ( kind = 4 ) l1

  real ( kind = crm_rknd ) ,intent(in   ) :: cc(in1,ido,4,l1)
  real ( kind = crm_rknd ) ,intent(inout) :: ch(in2,ido,l1,4)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) idp2
  integer ( kind = 4 ) im1
  integer ( kind = 4 ) im2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) m1
  integer ( kind = 4 ) m1d
  integer ( kind = 4 ) m2
  integer ( kind = 4 ) m2s
  real ( kind = crm_rknd ) sqrt2
  real ( kind = crm_rknd ) wa1(ido)
  real ( kind = crm_rknd ) wa2(ido)
  real ( kind = crm_rknd ) wa3(ido)

  m1d = ( m - 1 ) * im1 + 1
  m2s = 1 - im2
  sqrt2 = sqrt ( 2.0E+00 )

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

  if ( ido < 2 ) then
    return
  end if

  if ( 2 < ido ) then

    idp2 = ido + 2

    do k = 1, l1
      do i = 3, ido, 2
        ic = idp2 - i
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ch(m2,i-1,k,1) = (cc(m1,i-1,1,k)+cc(m1,ic-1,4,k)) &
            +(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k))
          ch(m2,i,k,1) = (cc(m1,i,1,k)-cc(m1,ic,4,k)) &
            +(cc(m1,i,3,k)-cc(m1,ic,2,k))
          ch(m2,i-1,k,2) = wa1(i-2)*((cc(m1,i-1,1,k)-cc(m1,ic-1,4,k)) &
            -(cc(m1,i,3,k)+cc(m1,ic,2,k)))-wa1(i-1) &
            *((cc(m1,i,1,k)+cc(m1,ic,4,k))+(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k)))
          ch(m2,i,k,2) = wa1(i-2)*((cc(m1,i,1,k)+cc(m1,ic,4,k)) &
            +(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))) + wa1(i-1) &
            *((cc(m1,i-1,1,k)-cc(m1,ic-1,4,k))-(cc(m1,i,3,k)+cc(m1,ic,2,k)))
          ch(m2,i-1,k,3) = wa2(i-2)*((cc(m1,i-1,1,k)+cc(m1,ic-1,4,k)) &
            -(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k))) - wa2(i-1) &
            *((cc(m1,i,1,k)-cc(m1,ic,4,k))-(cc(m1,i,3,k)-cc(m1,ic,2,k)))
          ch(m2,i,k,3) = wa2(i-2)*((cc(m1,i,1,k)-cc(m1,ic,4,k)) &
            -(cc(m1,i,3,k)-cc(m1,ic,2,k))) + wa2(i-1) &
            *((cc(m1,i-1,1,k)+cc(m1,ic-1,4,k))-(cc(m1,i-1,3,k) &
            +cc(m1,ic-1,2,k)))
          ch(m2,i-1,k,4) = wa3(i-2)*((cc(m1,i-1,1,k)-cc(m1,ic-1,4,k)) &
            +(cc(m1,i,3,k)+cc(m1,ic,2,k))) - wa3(i-1) &
            *((cc(m1,i,1,k)+cc(m1,ic,4,k))-(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k)))
          ch(m2,i,k,4) = wa3(i-2)*((cc(m1,i,1,k)+cc(m1,ic,4,k)) &
            -(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))) + wa3(i-1) &
            *((cc(m1,i-1,1,k)-cc(m1,ic-1,4,k))+(cc(m1,i,3,k)+cc(m1,ic,2,k)))
        end do
      end do
    end do

    if ( mod ( ido, 2 ) == 1 ) then
      return
    end if

  end if

  do k = 1, l1
    m2 = m2s
    do m1 = 1, m1d, im1
      m2 = m2 + im2
      ch(m2,ido,k,1) = (cc(m1,ido,1,k)+cc(m1,ido,3,k)) &
        +(cc(m1,ido,1,k)+cc(m1,ido,3,k))
      ch(m2,ido,k,2) = sqrt2*((cc(m1,ido,1,k)-cc(m1,ido,3,k)) &
        -(cc(m1,1,2,k)+cc(m1,1,4,k)))
      ch(m2,ido,k,3) = (cc(m1,1,4,k)-cc(m1,1,2,k)) &
        +(cc(m1,1,4,k)-cc(m1,1,2,k))
      ch(m2,ido,k,4) = -sqrt2*((cc(m1,ido,1,k)-cc(m1,ido,3,k)) &
        +(cc(m1,1,2,k)+cc(m1,1,4,k)))
    end do
  end do

  return
end
!========================================================================================
! MRADB5 is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mradb5 (m,ido,l1,cc,im1,in1,ch,im2,in2,wa1,wa2,wa3,wa4)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) ido
  integer ( kind = 4 ) in1
  integer ( kind = 4 ) in2
  integer ( kind = 4 ) l1

  real ( kind = crm_rknd ) arg
  real ( kind = crm_rknd ) ,intent(in   ) :: cc(in1,ido,5,l1)
  real ( kind = crm_rknd ) ,intent(inout) :: ch(in2,ido,l1,5)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) idp2
  integer ( kind = 4 ) im1
  integer ( kind = 4 ) im2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) m1
  integer ( kind = 4 ) m1d
  integer ( kind = 4 ) m2
  integer ( kind = 4 ) m2s
  real ( kind = crm_rknd ) ti11
  real ( kind = crm_rknd ) ti12
  real ( kind = crm_rknd ) tr11
  real ( kind = crm_rknd ) tr12
  real ( kind = crm_rknd ) wa1(ido)
  real ( kind = crm_rknd ) wa2(ido)
  real ( kind = crm_rknd ) wa3(ido)
  real ( kind = crm_rknd ) wa4(ido)

  m1d = ( m - 1 ) * im1 + 1
  m2s = 1 - im2
  arg = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 ) / 5.E+00
  tr11 = cos ( arg )
  ti11 = sin ( arg )
  tr12 = cos ( 2.0E+00 * arg )
  ti12 = sin ( 2.0E+00 * arg )

  do k = 1, l1
    m2 = m2s
    do m1 = 1, m1d, im1
      m2 = m2 + im2
      ch(m2,1,k,1) = cc(m1,1,1,k) + 2.0E+00 * cc(m1,ido,2,k) &
        + 2.0E+00 * cc(m1,ido,4,k)
      ch(m2,1,k,2) = ( cc(m1,1,1,k) + tr11 * 2.0E+00 * cc(m1,ido,2,k) &
        + tr12 * 2.0E+00 * cc(m1,ido,4,k) ) - ( ti11 * 2.0E+00 * cc(m1,1,3,k) &
        + ti12 * 2.0E+00 * cc(m1,1,5,k) )
      ch(m2,1,k,3) = ( cc(m1,1,1,k) + tr12 * 2.0E+00 * cc(m1,ido,2,k) &
        + tr11 * 2.0E+00 * cc(m1,ido,4,k) ) - ( ti12 * 2.0E+00 * cc(m1,1,3,k) &
        - ti11 * 2.0E+00 * cc(m1,1,5,k) )
      ch(m2,1,k,4) = ( cc(m1,1,1,k) + tr12 * 2.0E+00 * cc(m1,ido,2,k) &
        + tr11 * 2.0E+00 * cc(m1,ido,4,k) ) + ( ti12 * 2.0E+00 * cc(m1,1,3,k) &
        - ti11 * 2.0E+00 * cc(m1,1,5,k) )
      ch(m2,1,k,5) = ( cc(m1,1,1,k) + tr11 * 2.0E+00 * cc(m1,ido,2,k) &
        + tr12 * 2.0E+00 * cc(m1,ido,4,k) ) + ( ti11 * 2.0E+00 * cc(m1,1,3,k) &
        + ti12 * 2.0E+00 * cc(m1,1,5,k) )
    end do
  end do

  if ( ido == 1 ) then
    return
  end if

  idp2 = ido + 2
  do k = 1, l1
    do i = 3, ido, 2
      ic = idp2 - i
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch(m2,i-1,k,1) = cc(m1,i-1,1,k)+(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
          +(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k))
        ch(m2,i,k,1) = cc(m1,i,1,k)+(cc(m1,i,3,k)-cc(m1,ic,2,k)) &
          +(cc(m1,i,5,k)-cc(m1,ic,4,k))
        ch(m2,i-1,k,2) = wa1(i-2)*((cc(m1,i-1,1,k)+tr11* &
          (cc(m1,i-1,3,k)+cc(m1,ic-1,2,k))+tr12 &
          *(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))-(ti11*(cc(m1,i,3,k) &
          +cc(m1,ic,2,k))+ti12*(cc(m1,i,5,k)+cc(m1,ic,4,k)))) &
          -wa1(i-1)*((cc(m1,i,1,k)+tr11*(cc(m1,i,3,k)-cc(m1,ic,2,k)) &
          +tr12*(cc(m1,i,5,k)-cc(m1,ic,4,k)))+(ti11*(cc(m1,i-1,3,k) &
          -cc(m1,ic-1,2,k))+ti12*(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k))))
        ch(m2,i,k,2) = wa1(i-2)*((cc(m1,i,1,k)+tr11*(cc(m1,i,3,k) &
          -cc(m1,ic,2,k))+tr12*(cc(m1,i,5,k)-cc(m1,ic,4,k))) &
          +(ti11*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))+ti12 &
          *(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k)))) + wa1(i-1) &
          *((cc(m1,i-1,1,k)+tr11*(cc(m1,i-1,3,k) &
          +cc(m1,ic-1,2,k))+tr12*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k))) &
          -(ti11*(cc(m1,i,3,k)+cc(m1,ic,2,k))+ti12 &
          *(cc(m1,i,5,k)+cc(m1,ic,4,k))))
        ch(m2,i-1,k,3) = wa2(i-2) &
          *((cc(m1,i-1,1,k)+tr12*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
          +tr11*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))-(ti12*(cc(m1,i,3,k) &
          +cc(m1,ic,2,k))-ti11*(cc(m1,i,5,k)+cc(m1,ic,4,k)))) &
          -wa2(i-1) &
          *((cc(m1,i,1,k)+tr12*(cc(m1,i,3,k)- &
        cc(m1,ic,2,k))+tr11*(cc(m1,i,5,k)-cc(m1,ic,4,k))) &
          +(ti12*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))-ti11 &
          *(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k))))
        ch(m2,i,k,3) = wa2(i-2) &
          *((cc(m1,i,1,k)+tr12*(cc(m1,i,3,k)- &
          cc(m1,ic,2,k))+tr11*(cc(m1,i,5,k)-cc(m1,ic,4,k))) &
          +(ti12*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))-ti11 &
          *(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k)))) &
          + wa2(i-1) &
          *((cc(m1,i-1,1,k)+tr12*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
          +tr11*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))-(ti12*(cc(m1,i,3,k) &
          +cc(m1,ic,2,k))-ti11*(cc(m1,i,5,k)+cc(m1,ic,4,k))))
        ch(m2,i-1,k,4) = wa3(i-2) &
          *((cc(m1,i-1,1,k)+tr12*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
          +tr11*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))+(ti12*(cc(m1,i,3,k) &
          +cc(m1,ic,2,k))-ti11*(cc(m1,i,5,k)+cc(m1,ic,4,k)))) &
          -wa3(i-1) &
          *((cc(m1,i,1,k)+tr12*(cc(m1,i,3,k)- &
          cc(m1,ic,2,k))+tr11*(cc(m1,i,5,k)-cc(m1,ic,4,k))) &
          -(ti12*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))-ti11 &
          *(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k))))
        ch(m2,i,k,4) = wa3(i-2) &
          *((cc(m1,i,1,k)+tr12*(cc(m1,i,3,k)- &
          cc(m1,ic,2,k))+tr11*(cc(m1,i,5,k)-cc(m1,ic,4,k))) &
          -(ti12*(cc(m1,i-1,3,k)-cc(m1,ic-1,2,k))-ti11 &
          *(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k)))) &
          + wa3(i-1) &
          *((cc(m1,i-1,1,k)+tr12*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
          +tr11*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))+(ti12*(cc(m1,i,3,k) &
          +cc(m1,ic,2,k))-ti11*(cc(m1,i,5,k)+cc(m1,ic,4,k))))
        ch(m2,i-1,k,5) = wa4(i-2) &
          *((cc(m1,i-1,1,k)+tr11*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
          +tr12*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))+(ti11*(cc(m1,i,3,k) &
          +cc(m1,ic,2,k))+ti12*(cc(m1,i,5,k)+cc(m1,ic,4,k)))) &
          -wa4(i-1) &
          *((cc(m1,i,1,k)+tr11*(cc(m1,i,3,k)-cc(m1,ic,2,k)) &
          +tr12*(cc(m1,i,5,k)-cc(m1,ic,4,k)))-(ti11*(cc(m1,i-1,3,k) &
          -cc(m1,ic-1,2,k))+ti12*(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k))))
        ch(m2,i,k,5) = wa4(i-2) &
          *((cc(m1,i,1,k)+tr11*(cc(m1,i,3,k)-cc(m1,ic,2,k)) &
          +tr12*(cc(m1,i,5,k)-cc(m1,ic,4,k)))-(ti11*(cc(m1,i-1,3,k) &
          -cc(m1,ic-1,2,k))+ti12*(cc(m1,i-1,5,k)-cc(m1,ic-1,4,k)))) &
          + wa4(i-1) &
          *((cc(m1,i-1,1,k)+tr11*(cc(m1,i-1,3,k)+cc(m1,ic-1,2,k)) &
          +tr12*(cc(m1,i-1,5,k)+cc(m1,ic-1,4,k)))+(ti11*(cc(m1,i,3,k) &
          +cc(m1,ic,2,k))+ti12*(cc(m1,i,5,k)+cc(m1,ic,4,k))))
      end do
    end do
  end do

  return
end
!========================================================================================
! MRADBG is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mradbg (m,ido,ip,l1,idl1,cc,c1,c2,im1,in1,ch,ch2,im2,in2,wa)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) idl1
  integer ( kind = 4 ) ido
  integer ( kind = 4 ) in1
  integer ( kind = 4 ) in2
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) l1

  real ( kind = crm_rknd ) ai1
  real ( kind = crm_rknd ) ai2
  real ( kind = crm_rknd ) ar1
  real ( kind = crm_rknd ) ar1h
  real ( kind = crm_rknd ) ar2
  real ( kind = crm_rknd ) ar2h
  real ( kind = crm_rknd ) arg
  real ( kind = crm_rknd ) c1(in1,ido,l1,ip)
  real ( kind = crm_rknd ) c2(in1,idl1,ip)
  real ( kind = crm_rknd ) cc(in1,ido,ip,l1) 
  real ( kind = crm_rknd ) ,intent(inout) :: ch(in2,ido,l1,ip)
  real ( kind = crm_rknd ) ch2(in2,idl1,ip)
  real ( kind = crm_rknd ) dc2
  real ( kind = crm_rknd ) dcp
  real ( kind = crm_rknd ) ds2
  real ( kind = crm_rknd ) dsp
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) idij
  integer ( kind = 4 ) idp2
  integer ( kind = 4 ) ik
  integer ( kind = 4 ) im1
  integer ( kind = 4 ) im2
  integer ( kind = 4 ) ipp2
  integer ( kind = 4 ) ipph
  integer ( kind = 4 ) is
  integer ( kind = 4 ) j
  integer ( kind = 4 ) j2
  integer ( kind = 4 ) jc
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  integer ( kind = 4 ) lc
  integer ( kind = 4 ) m
  integer ( kind = 4 ) m1
  integer ( kind = 4 ) m1d
  integer ( kind = 4 ) m2
  integer ( kind = 4 ) m2s
  integer ( kind = 4 ) nbd
  real ( kind = crm_rknd ) tpi
  real ( kind = crm_rknd ) wa(ido)

  m1d = ( m - 1 ) * im1 + 1
  m2s = 1 - im2
  tpi = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 )
  arg = tpi / real ( ip, kind = crm_rknd )
  dcp = cos ( arg )
  dsp = sin ( arg )
  idp2 = ido + 2
  nbd = ( ido - 1 ) / 2
  ipp2 = ip + 2
  ipph = ( ip + 1 ) / 2

  if ( ido < l1 ) then

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

  else

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

  end if

  do j = 2, ipph
    jc = ipp2 - j
    j2 = j + j
    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch(m2,1,k,j) = cc(m1,ido,j2-2,k) + cc(m1,ido,j2-2,k)
        ch(m2,1,k,jc) = cc(m1,1,j2-1,k) + cc(m1,1,j2-1,k)
      end do
    end do
  end do

  if ( ido == 1 ) then

  else if ( nbd < l1 ) then

    do j = 2, ipph
      jc = ipp2 - j
      do i = 3, ido, 2
        ic = idp2 - i
        do k = 1, l1
          m2 = m2s
          do m1 = 1, m1d, im1
            m2 = m2 + im2
            ch(m2,i-1,k,j)  = cc(m1,i-1,2*j-1,k) + cc(m1,ic-1,2*j-2,k)
            ch(m2,i-1,k,jc) = cc(m1,i-1,2*j-1,k) - cc(m1,ic-1,2*j-2,k)
            ch(m2,i,k,j)    = cc(m1,i,2*j-1,k)   - cc(m1,ic,2*j-2,k)
            ch(m2,i,k,jc)   = cc(m1,i,2*j-1,k)   + cc(m1,ic,2*j-2,k)
          end do
        end do
      end do
    end do

  else

    do j = 2, ipph
      jc = ipp2 - j
      do k = 1, l1
        do i = 3, ido, 2
          ic = idp2 - i
          m2 = m2s
          do m1 = 1, m1d, im1
            m2 = m2 + im2
            ch(m2,i-1,k,j)  = cc(m1,i-1,2*j-1,k) + cc(m1,ic-1,2*j-2,k)
            ch(m2,i-1,k,jc) = cc(m1,i-1,2*j-1,k) - cc(m1,ic-1,2*j-2,k)
            ch(m2,i,k,j)    = cc(m1,i,2*j-1,k)   - cc(m1,ic,2*j-2,k)
            ch(m2,i,k,jc)   = cc(m1,i,2*j-1,k)   + cc(m1,ic,2*j-2,k)
          end do
        end do
      end do
    end do

  end if

  ar1 = 1.0E+00
  ai1 = 0.0E+00

  do l = 2, ipph

    lc = ipp2 - l
    ar1h = dcp * ar1 - dsp * ai1
    ai1 =  dcp * ai1 + dsp * ar1
    ar1 = ar1h

    do ik = 1, idl1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        c2(m1,ik,l)  = ch2(m2,ik,1) + ar1 * ch2(m2,ik,2)
        c2(m1,ik,lc) =                ai1 * ch2(m2,ik,ip)
      end do
    end do

    dc2 = ar1
    ds2 = ai1
    ar2 = ar1
    ai2 = ai1

    do j = 3, ipph
      jc = ipp2 - j
      ar2h = dc2 * ar2 - ds2 * ai2
      ai2  = dc2 * ai2 + ds2 * ar2
      ar2 = ar2h
      do ik = 1, idl1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          c2(m1,ik,l)  = c2(m1,ik,l)  + ar2 * ch2(m2,ik,j)
          c2(m1,ik,lc) = c2(m1,ik,lc) + ai2 * ch2(m2,ik,jc)
        end do
      end do
    end do

  end do

  do j = 2, ipph
    do ik = 1, idl1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch2(m2,ik,1) = ch2(m2,ik,1) + ch2(m2,ik,j)
      end do
    end do
  end do

  do j = 2, ipph
    jc = ipp2 - j
    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch(m2,1,k,j)  = c1(m1,1,k,j) - c1(m1,1,k,jc)
        ch(m2,1,k,jc) = c1(m1,1,k,j) + c1(m1,1,k,jc)
      end do
    end do
  end do

  if ( ido == 1 ) then

  else if ( nbd < l1 ) then

    do j = 2, ipph
      jc = ipp2 - j
      do i = 3, ido, 2
        do k = 1, l1
          m2 = m2s
          do m1 = 1, m1d, im1
            m2 = m2 + im2
            ch(m2,i-1,k,j)  = c1(m1,i-1,k,j) - c1(m1,i,k,jc)
            ch(m2,i-1,k,jc) = c1(m1,i-1,k,j) + c1(m1,i,k,jc)
            ch(m2,i,k,j)    = c1(m1,i,k,j)   + c1(m1,i-1,k,jc)
            ch(m2,i,k,jc)   = c1(m1,i,k,j)   - c1(m1,i-1,k,jc)
          end do
        end do
      end do
    end do

  else

    do j = 2, ipph
      jc = ipp2 - j
      do k = 1, l1
        do i = 3, ido, 2
          m2 = m2s
          do m1 = 1, m1d, im1
            m2 = m2 + im2
            ch(m2,i-1,k,j)  = c1(m1,i-1,k,j) - c1(m1,i,k,jc)
            ch(m2,i-1,k,jc) = c1(m1,i-1,k,j) + c1(m1,i,k,jc)
            ch(m2,i,k,j)    = c1(m1,i,k,j)   + c1(m1,i-1,k,jc)
            ch(m2,i,k,jc)   = c1(m1,i,k,j)   - c1(m1,i-1,k,jc)
          end do
        end do
      end do
    end do

  end if

  if ( ido == 1 ) then
    return
  end if

  do ik = 1, idl1
    m2 = m2s
    do m1 = 1, m1d, im1
      m2 = m2 + im2
      c2(m1,ik,1) = ch2(m2,ik,1)
    end do
  end do

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

  if ( l1 < nbd ) then

    is = -ido

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

  else

    is = -ido

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

  end if

  return
end
!========================================================================================
! MRADF2 is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mradf2 (m,ido,l1,cc,im1,in1,ch,im2,in2,wa1)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) ido
  integer ( kind = 4 ) in1
  integer ( kind = 4 ) in2
  integer ( kind = 4 ) l1

  real ( kind = crm_rknd ) cc(in1,ido,l1,2)
  real ( kind = crm_rknd ) ,intent(inout) :: ch(in2,ido,2,l1)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) idp2
  integer ( kind = 4 ) im1
  integer ( kind = 4 ) im2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) m1
  integer ( kind = 4 ) m1d
  integer ( kind = 4 ) m2
  integer ( kind = 4 ) m2s
  real ( kind = crm_rknd ) wa1(ido)

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

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

  if ( ido < 2 ) then
    return
  end if

  if ( 2 < ido ) then

    idp2 = ido + 2

    do k = 1, l1
      do i = 3, ido, 2
        ic = idp2 - i
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ch(m2,i,1,k) =    cc(m1,i,k,1)   + ( wa1(i-2) * cc(m1,i,k,2) &
                                             - wa1(i-1) * cc(m1,i-1,k,2) )
          ch(m2,ic,2,k)  = -cc(m1,i,k,1)   + ( wa1(i-2) * cc(m1,i,k,2) &
                                             - wa1(i-1) * cc(m1,i-1,k,2) )
          ch(m2,i-1,1,k)  = cc(m1,i-1,k,1) + ( wa1(i-2) * cc(m1,i-1,k,2) &
                                             + wa1(i-1) * cc(m1,i,k,2))
          ch(m2,ic-1,2,k) = cc(m1,i-1,k,1) - ( wa1(i-2) * cc(m1,i-1,k,2) &
                                             + wa1(i-1) * cc(m1,i,k,2) )
        end do
      end do
    end do

    if ( mod ( ido, 2 ) == 1 ) then
      return
    end if

  end if

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

  return
end
!========================================================================================
! MRADF3 is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mradf3 (m,ido,l1,cc,im1,in1,ch,im2,in2,wa1,wa2)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) ido
  integer ( kind = 4 ) in1
  integer ( kind = 4 ) in2
  integer ( kind = 4 ) l1

  real ( kind = crm_rknd ) arg
  real ( kind = crm_rknd ) cc(in1,ido,l1,3)
  real ( kind = crm_rknd ) ,intent(inout) :: ch(in2,ido,3,l1)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) idp2
  integer ( kind = 4 ) im1
  integer ( kind = 4 ) im2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) m1
  integer ( kind = 4 ) m1d
  integer ( kind = 4 ) m2
  integer ( kind = 4 ) m2s
  real ( kind = crm_rknd ) taui
  real ( kind = crm_rknd ) taur
  real ( kind = crm_rknd ) wa1(ido)
  real ( kind = crm_rknd ) wa2(ido)

  m1d = ( m - 1 ) * im1 + 1
  m2s = 1 - im2
  arg = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 ) / 3.0E+00
  taur = cos ( arg )
  taui = sin ( arg )

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

  if ( ido == 1 ) then
    return
  end if

  idp2 = ido + 2

  do k = 1, l1
    do i = 3, ido, 2
      ic = idp2 - i
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch(m2,i-1,1,k) = cc(m1,i-1,k,1)+((wa1(i-2)*cc(m1,i-1,k,2)+ &
          wa1(i-1)*cc(m1,i,k,2))+(wa2(i-2)*cc(m1,i-1,k,3) + wa2(i-1)* &
          cc(m1,i,k,3)))
        ch(m2,i,1,k) = cc(m1,i,k,1)+((wa1(i-2)*cc(m1,i,k,2)- &
          wa1(i-1)*cc(m1,i-1,k,2))+(wa2(i-2)*cc(m1,i,k,3)-wa2(i-1)* &
          cc(m1,i-1,k,3)))
        ch(m2,i-1,3,k) = (cc(m1,i-1,k,1)+taur*((wa1(i-2)* &
          cc(m1,i-1,k,2) + wa1(i-1)*cc(m1,i,k,2))+(wa2(i-2)* &
          cc(m1,i-1,k,3) + wa2(i-1)*cc(m1,i,k,3))))+(taui*((wa1(i-2)* &
          cc(m1,i,k,2) - wa1(i-1)*cc(m1,i-1,k,2))-(wa2(i-2)* &
          cc(m1,i,k,3) - wa2(i-1)*cc(m1,i-1,k,3))))
        ch(m2,ic-1,2,k) = (cc(m1,i-1,k,1)+taur*((wa1(i-2)* &
          cc(m1,i-1,k,2) + wa1(i-1)*cc(m1,i,k,2))+(wa2(i-2)* &
          cc(m1,i-1,k,3) + wa2(i-1)*cc(m1,i,k,3))))-(taui*((wa1(i-2)* &
          cc(m1,i,k,2) - wa1(i-1)*cc(m1,i-1,k,2)) - ( wa2(i-2)* &
          cc(m1,i,k,3)-wa2(i-1)*cc(m1,i-1,k,3))))
        ch(m2,i,3,k) = (cc(m1,i,k,1)+taur*((wa1(i-2)*cc(m1,i,k,2)- &
          wa1(i-1)*cc(m1,i-1,k,2))+(wa2(i-2)*cc(m1,i,k,3) - wa2(i-1)* &
          cc(m1,i-1,k,3))))+(taui*((wa2(i-2)*cc(m1,i-1,k,3) + wa2(i-1)* &
          cc(m1,i,k,3))-(wa1(i-2)*cc(m1,i-1,k,2) + wa1(i-1)* &
          cc(m1,i,k,2))))
        ch(m2,ic,2,k) = (taui*((wa2(i-2)*cc(m1,i-1,k,3) + wa2(i-1)* &
          cc(m1,i,k,3))-(wa1(i-2)*cc(m1,i-1,k,2) + wa1(i-1)* &
          cc(m1,i,k,2))))-(cc(m1,i,k,1)+taur*((wa1(i-2)*cc(m1,i,k,2)- &
          wa1(i-1)*cc(m1,i-1,k,2))+(wa2(i-2)*cc(m1,i,k,3)-wa2(i-1)* &
          cc(m1,i-1,k,3))))
      end do
    end do
  end do

  return
end
!========================================================================================
! MRADF4 is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mradf4 (m,ido,l1,cc,im1,in1,ch,im2,in2,wa1,wa2,wa3)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) ido
  integer ( kind = 4 ) in1
  integer ( kind = 4 ) in2
  integer ( kind = 4 ) l1

  real    ( kind = crm_rknd ) ,intent(in   ) :: cc(in1,ido,l1,4)
  real    ( kind = crm_rknd ) ,intent(inout) :: ch(in2,ido,4,l1)
  real    ( kind = crm_rknd ) ,intent(inout) :: wa1(ido)
  real    ( kind = crm_rknd ) ,intent(inout) :: wa2(ido)
  real    ( kind = crm_rknd ) ,intent(inout) :: wa3(ido)
  real    ( kind = crm_rknd ) hsqt2
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) idp2
  integer ( kind = 4 ) im1
  integer ( kind = 4 ) im2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) m1
  integer ( kind = 4 ) m1d
  integer ( kind = 4 ) m2
  integer ( kind = 4 ) m2s
  

! m   = lot  = nz
! im1 = (jump = nx) or (1)
! im2 = (1)         or (jump = nx)
! in1 = (inc = 1)   or (lot = nz)
! in2 = (lot = nz)  or (inc = 1)

! m1d = { (nz-1)*nx+1 } or { (nz-1)*1+1 = nz }

  hsqt2 = sqrt ( 2.0E+00 ) / 2.0E+00
  m1d = ( m - 1 ) * im1 + 1
  m2s = 1 - im2

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

  if ( ido < 2 ) then
    return
  end if

  if ( 2 < ido ) then

    idp2 = ido + 2

    do k = 1, l1
      do i = 3, ido, 2
        ic = idp2 - i
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ch(m2,i-1,1,k) = ((wa1(i-2)*cc(m1,i-1,k,2) + wa1(i-1)* &
            cc(m1,i,k,2))+(wa3(i-2)*cc(m1,i-1,k,4) + wa3(i-1)* &
            cc(m1,i,k,4)))+(cc(m1,i-1,k,1)+(wa2(i-2)*cc(m1,i-1,k,3)+ &
            wa2(i-1)*cc(m1,i,k,3)))
          ch(m2,ic-1,4,k) = (cc(m1,i-1,k,1)+(wa2(i-2)*cc(m1,i-1,k,3)+ &
            wa2(i-1)*cc(m1,i,k,3)))-((wa1(i-2)*cc(m1,i-1,k,2)+ &
            wa1(i-1)*cc(m1,i,k,2))+(wa3(i-2)*cc(m1,i-1,k,4)+ &
            wa3(i-1)*cc(m1,i,k,4)))
          ch(m2,i,1,k) = ((wa1(i-2)*cc(m1,i,k,2)-wa1(i-1)* &
            cc(m1,i-1,k,2))+(wa3(i-2)*cc(m1,i,k,4)-wa3(i-1)* &
            cc(m1,i-1,k,4)))+(cc(m1,i,k,1)+(wa2(i-2)*cc(m1,i,k,3)- &
            wa2(i-1)*cc(m1,i-1,k,3)))
          ch(m2,ic,4,k) = ((wa1(i-2)*cc(m1,i,k,2)-wa1(i-1)* &
            cc(m1,i-1,k,2))+(wa3(i-2)*cc(m1,i,k,4)-wa3(i-1)* &
            cc(m1,i-1,k,4)))-(cc(m1,i,k,1)+(wa2(i-2)*cc(m1,i,k,3)- &
            wa2(i-1)*cc(m1,i-1,k,3)))
          ch(m2,i-1,3,k) = ((wa1(i-2)*cc(m1,i,k,2)-wa1(i-1)* &
            cc(m1,i-1,k,2))-(wa3(i-2)*cc(m1,i,k,4)-wa3(i-1)* &
            cc(m1,i-1,k,4)))+(cc(m1,i-1,k,1)-(wa2(i-2)*cc(m1,i-1,k,3)+ &
            wa2(i-1)*cc(m1,i,k,3)))
          ch(m2,ic-1,2,k) = (cc(m1,i-1,k,1)-(wa2(i-2)*cc(m1,i-1,k,3)+ &
            wa2(i-1)*cc(m1,i,k,3)))-((wa1(i-2)*cc(m1,i,k,2)-wa1(i-1)* &
            cc(m1,i-1,k,2))-(wa3(i-2)*cc(m1,i,k,4)-wa3(i-1)* &
            cc(m1,i-1,k,4)))
          ch(m2,i,3,k) = ((wa3(i-2)*cc(m1,i-1,k,4) + wa3(i-1)* &
            cc(m1,i,k,4))-(wa1(i-2)*cc(m1,i-1,k,2) + wa1(i-1)* &
            cc(m1,i,k,2)))+(cc(m1,i,k,1)-(wa2(i-2)*cc(m1,i,k,3)- &
            wa2(i-1)*cc(m1,i-1,k,3)))
          ch(m2,ic,2,k) = ((wa3(i-2)*cc(m1,i-1,k,4) + wa3(i-1)* &
            cc(m1,i,k,4))-(wa1(i-2)*cc(m1,i-1,k,2) + wa1(i-1)* &
            cc(m1,i,k,2)))-(cc(m1,i,k,1)-(wa2(i-2)*cc(m1,i,k,3)- &
            wa2(i-1)*cc(m1,i-1,k,3)))
        end do
      end do
    end do

    if ( mod ( ido, 2 ) == 1 ) then
      return
    end if

  end if

  do k = 1, l1
    m2 = m2s
    do m1 = 1, m1d, im1
      m2 = m2 + im2
      ch(m2,ido,1,k) = cc(m1,ido,k,1) &
        + (  hsqt2 * ( cc(m1,ido,k,2) - cc(m1,ido,k,4) ) )
      ch(m2,ido,3,k) = cc(m1,ido,k,1) &
        - (  hsqt2 * ( cc(m1,ido,k,2) - cc(m1,ido,k,4) ) )
      ch(m2,1,2,k) =  -cc(m1,ido,k,3) &
        + ( -hsqt2 * ( cc(m1,ido,k,2) + cc(m1,ido,k,4) ) )
      ch(m2,1,4,k) =   cc(m1,ido,k,3) &
        + ( -hsqt2 * ( cc(m1,ido,k,2) + cc(m1,ido,k,4) ) )
    end do
  end do

  return
end
!========================================================================================
! MRADF5 is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mradf5 (m,ido,l1,cc,im1,in1,ch,im2,in2,wa1,wa2,wa3,wa4)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) ido
  integer ( kind = 4 ) in1
  integer ( kind = 4 ) in2
  integer ( kind = 4 ) l1

  real ( kind = crm_rknd ) arg
  real ( kind = crm_rknd ) cc(in1,ido,l1,5)
  real ( kind = crm_rknd ) ,intent(inout) :: ch(in2,ido,5,l1)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) idp2
  integer ( kind = 4 ) im1
  integer ( kind = 4 ) im2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) m1
  integer ( kind = 4 ) m1d
  integer ( kind = 4 ) m2
  integer ( kind = 4 ) m2s
  real ( kind = crm_rknd ) ti11
  real ( kind = crm_rknd ) ti12
  real ( kind = crm_rknd ) tr11
  real ( kind = crm_rknd ) tr12
  real ( kind = crm_rknd ) wa1(ido)
  real ( kind = crm_rknd ) wa2(ido)
  real ( kind = crm_rknd ) wa3(ido)
  real ( kind = crm_rknd ) wa4(ido)

  m1d = ( m - 1 ) * im1 + 1
  m2s = 1 - im2
  arg = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 ) / 5.0E+00
  tr11 = cos ( arg )
  ti11 = sin ( arg )
  tr12 = cos ( 2.0E+00 * arg )
  ti12 = sin ( 2.0E+00 * arg )

  do k = 1, l1
    m2 = m2s
    do m1 = 1, m1d, im1
      m2 = m2 + im2
      ch(m2,1,1,k) = cc(m1,1,k,1) + ( cc(m1,1,k,5) + cc(m1,1,k,2) ) &
                                  + ( cc(m1,1,k,4) + cc(m1,1,k,3) )
      ch(m2,ido,2,k) = cc(m1,1,k,1) + tr11 * ( cc(m1,1,k,5) + cc(m1,1,k,2) ) &
                                    + tr12 * ( cc(m1,1,k,4) + cc(m1,1,k,3) )
      ch(m2,1,3,k) = ti11 * ( cc(m1,1,k,5) - cc(m1,1,k,2) ) &
                   + ti12 * ( cc(m1,1,k,4) - cc(m1,1,k,3) )
      ch(m2,ido,4,k) = cc(m1,1,k,1) + tr12 * ( cc(m1,1,k,5) + cc(m1,1,k,2) ) &
                                    + tr11 * ( cc(m1,1,k,4) + cc(m1,1,k,3) )
      ch(m2,1,5,k) = ti12 * ( cc(m1,1,k,5) - cc(m1,1,k,2) ) &
                   - ti11 * ( cc(m1,1,k,4) - cc(m1,1,k,3) )
    end do
  end do

  if ( ido == 1 ) then
    return
  end if

  idp2 = ido + 2

  do k = 1, l1
    do i = 3, ido, 2
      ic = idp2 - i
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch(m2,i-1,1,k) = cc(m1,i-1,k,1)+((wa1(i-2)*cc(m1,i-1,k,2)+ &
          wa1(i-1)*cc(m1,i,k,2))+(wa4(i-2)*cc(m1,i-1,k,5)+wa4(i-1)* &
          cc(m1,i,k,5)))+((wa2(i-2)*cc(m1,i-1,k,3)+wa2(i-1)* &
          cc(m1,i,k,3))+(wa3(i-2)*cc(m1,i-1,k,4)+ &
          wa3(i-1)*cc(m1,i,k,4)))
        ch(m2,i,1,k) = cc(m1,i,k,1)+((wa1(i-2)*cc(m1,i,k,2)- &
          wa1(i-1)*cc(m1,i-1,k,2))+(wa4(i-2)*cc(m1,i,k,5)-wa4(i-1)* &
          cc(m1,i-1,k,5)))+((wa2(i-2)*cc(m1,i,k,3)-wa2(i-1)* &
          cc(m1,i-1,k,3))+(wa3(i-2)*cc(m1,i,k,4)-wa3(i-1)* &
          cc(m1,i-1,k,4)))
        ch(m2,i-1,3,k) = cc(m1,i-1,k,1)+tr11* &
          ( wa1(i-2)*cc(m1,i-1,k,2)+wa1(i-1)*cc(m1,i,k,2) &
          +wa4(i-2)*cc(m1,i-1,k,5)+wa4(i-1)*cc(m1,i,k,5))+tr12* &
          ( wa2(i-2)*cc(m1,i-1,k,3)+wa2(i-1)*cc(m1,i,k,3) &
          +wa3(i-2)*cc(m1,i-1,k,4)+wa3(i-1)*cc(m1,i,k,4))+ti11* &
          ( wa1(i-2)*cc(m1,i,k,2)-wa1(i-1)*cc(m1,i-1,k,2) &
          -(wa4(i-2)*cc(m1,i,k,5)-wa4(i-1)*cc(m1,i-1,k,5)))+ti12* &
          ( wa2(i-2)*cc(m1,i,k,3)-wa2(i-1)*cc(m1,i-1,k,3) &
          -(wa3(i-2)*cc(m1,i,k,4)-wa3(i-1)*cc(m1,i-1,k,4)))
        ch(m2,ic-1,2,k) = cc(m1,i-1,k,1)+tr11* &
          ( wa1(i-2)*cc(m1,i-1,k,2)+wa1(i-1)*cc(m1,i,k,2) &
          +wa4(i-2)*cc(m1,i-1,k,5)+wa4(i-1)*cc(m1,i,k,5))+tr12* &
          ( wa2(i-2)*cc(m1,i-1,k,3)+wa2(i-1)*cc(m1,i,k,3) &
          +wa3(i-2)*cc(m1,i-1,k,4)+wa3(i-1)*cc(m1,i,k,4))-(ti11* &
          ( wa1(i-2)*cc(m1,i,k,2)-wa1(i-1)*cc(m1,i-1,k,2) &
          -(wa4(i-2)*cc(m1,i,k,5)-wa4(i-1)*cc(m1,i-1,k,5)))+ti12* &
          ( wa2(i-2)*cc(m1,i,k,3)-wa2(i-1)*cc(m1,i-1,k,3) &
          -(wa3(i-2)*cc(m1,i,k,4)-wa3(i-1)*cc(m1,i-1,k,4))))
        ch(m2,i,3,k) = (cc(m1,i,k,1)+tr11*((wa1(i-2)*cc(m1,i,k,2)- &
          wa1(i-1)*cc(m1,i-1,k,2))+(wa4(i-2)*cc(m1,i,k,5)-wa4(i-1)* &
          cc(m1,i-1,k,5)))+tr12*((wa2(i-2)*cc(m1,i,k,3)-wa2(i-1)* &
          cc(m1,i-1,k,3))+(wa3(i-2)*cc(m1,i,k,4)-wa3(i-1)* &
          cc(m1,i-1,k,4))))+(ti11*((wa4(i-2)*cc(m1,i-1,k,5)+ &
          wa4(i-1)*cc(m1,i,k,5))-(wa1(i-2)*cc(m1,i-1,k,2)+wa1(i-1)* &
          cc(m1,i,k,2)))+ti12*((wa3(i-2)*cc(m1,i-1,k,4)+wa3(i-1)* &
          cc(m1,i,k,4))-(wa2(i-2)*cc(m1,i-1,k,3)+wa2(i-1)* &
          cc(m1,i,k,3))))
        ch(m2,ic,2,k) = (ti11*((wa4(i-2)*cc(m1,i-1,k,5)+wa4(i-1)* &
          cc(m1,i,k,5))-(wa1(i-2)*cc(m1,i-1,k,2)+wa1(i-1)* &
          cc(m1,i,k,2)))+ti12*((wa3(i-2)*cc(m1,i-1,k,4)+wa3(i-1)* &
          cc(m1,i,k,4))-(wa2(i-2)*cc(m1,i-1,k,3)+wa2(i-1)* &
          cc(m1,i,k,3))))-(cc(m1,i,k,1)+tr11*((wa1(i-2)*cc(m1,i,k,2)- &
          wa1(i-1)*cc(m1,i-1,k,2))+(wa4(i-2)*cc(m1,i,k,5)-wa4(i-1)* &
          cc(m1,i-1,k,5)))+tr12*((wa2(i-2)*cc(m1,i,k,3)-wa2(i-1)* &
          cc(m1,i-1,k,3))+(wa3(i-2)*cc(m1,i,k,4)-wa3(i-1)* &
          cc(m1,i-1,k,4))))
        ch(m2,i-1,5,k) = (cc(m1,i-1,k,1)+tr12*((wa1(i-2)* &
          cc(m1,i-1,k,2)+wa1(i-1)*cc(m1,i,k,2))+(wa4(i-2)* &
          cc(m1,i-1,k,5)+wa4(i-1)*cc(m1,i,k,5)))+tr11*((wa2(i-2)* &
          cc(m1,i-1,k,3)+wa2(i-1)*cc(m1,i,k,3))+(wa3(i-2)* &
          cc(m1,i-1,k,4)+wa3(i-1)*cc(m1,i,k,4))))+(ti12*((wa1(i-2)* &
          cc(m1,i,k,2)-wa1(i-1)*cc(m1,i-1,k,2))-(wa4(i-2)* &
          cc(m1,i,k,5)-wa4(i-1)*cc(m1,i-1,k,5)))-ti11*((wa2(i-2)* &
          cc(m1,i,k,3)-wa2(i-1)*cc(m1,i-1,k,3))-(wa3(i-2)* &
          cc(m1,i,k,4)-wa3(i-1)*cc(m1,i-1,k,4))))
        ch(m2,ic-1,4,k) = (cc(m1,i-1,k,1)+tr12*((wa1(i-2)* &
          cc(m1,i-1,k,2)+wa1(i-1)*cc(m1,i,k,2))+(wa4(i-2)* &
          cc(m1,i-1,k,5)+wa4(i-1)*cc(m1,i,k,5)))+tr11*((wa2(i-2)* &
          cc(m1,i-1,k,3)+wa2(i-1)*cc(m1,i,k,3))+(wa3(i-2)* &
          cc(m1,i-1,k,4)+wa3(i-1)*cc(m1,i,k,4))))-(ti12*((wa1(i-2)* &
          cc(m1,i,k,2)-wa1(i-1)*cc(m1,i-1,k,2))-(wa4(i-2)* &
          cc(m1,i,k,5)-wa4(i-1)*cc(m1,i-1,k,5)))-ti11*((wa2(i-2)* &
          cc(m1,i,k,3)-wa2(i-1)*cc(m1,i-1,k,3))-(wa3(i-2)* &
          cc(m1,i,k,4)-wa3(i-1)*cc(m1,i-1,k,4))))
        ch(m2,i,5,k) = (cc(m1,i,k,1)+tr12*((wa1(i-2)*cc(m1,i,k,2)- &
          wa1(i-1)*cc(m1,i-1,k,2))+(wa4(i-2)*cc(m1,i,k,5)-wa4(i-1)* &
          cc(m1,i-1,k,5)))+tr11*((wa2(i-2)*cc(m1,i,k,3)-wa2(i-1)* &
          cc(m1,i-1,k,3))+(wa3(i-2)*cc(m1,i,k,4)-wa3(i-1)* &
          cc(m1,i-1,k,4))))+(ti12*((wa4(i-2)*cc(m1,i-1,k,5)+ &
          wa4(i-1)*cc(m1,i,k,5))-(wa1(i-2)*cc(m1,i-1,k,2)+wa1(i-1)* &
          cc(m1,i,k,2)))-ti11*((wa3(i-2)*cc(m1,i-1,k,4)+wa3(i-1)* &
          cc(m1,i,k,4))-(wa2(i-2)*cc(m1,i-1,k,3)+wa2(i-1)* &
          cc(m1,i,k,3))))
        ch(m2,ic,4,k) = (ti12*((wa4(i-2)*cc(m1,i-1,k,5)+wa4(i-1)* &
          cc(m1,i,k,5))-(wa1(i-2)*cc(m1,i-1,k,2)+wa1(i-1)* &
          cc(m1,i,k,2)))-ti11*((wa3(i-2)*cc(m1,i-1,k,4)+wa3(i-1)* &
          cc(m1,i,k,4))-(wa2(i-2)*cc(m1,i-1,k,3)+wa2(i-1)* &
          cc(m1,i,k,3))))-(cc(m1,i,k,1)+tr12*((wa1(i-2)*cc(m1,i,k,2)- &
          wa1(i-1)*cc(m1,i-1,k,2))+(wa4(i-2)*cc(m1,i,k,5)-wa4(i-1)* &
          cc(m1,i-1,k,5)))+tr11*((wa2(i-2)*cc(m1,i,k,3)-wa2(i-1)* &
          cc(m1,i-1,k,3))+(wa3(i-2)*cc(m1,i,k,4)-wa3(i-1)* &
          cc(m1,i-1,k,4))))
      end do
    end do
  end do

  return
end
!========================================================================================
! MRADFG is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mradfg (m,ido,ip,l1,idl1,cc,c1,c2,im1,in1,ch,ch2,im2,in2,wa)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) idl1
  integer ( kind = 4 ) ido
  integer ( kind = 4 ) in1
  integer ( kind = 4 ) in2
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) l1

  real ( kind = crm_rknd ) ai1
  real ( kind = crm_rknd ) ai2
  real ( kind = crm_rknd ) ar1
  real ( kind = crm_rknd ) ar1h
  real ( kind = crm_rknd ) ar2
  real ( kind = crm_rknd ) ar2h
  real ( kind = crm_rknd ) arg
  real ( kind = crm_rknd ) c1(in1,ido,l1,ip)
  real ( kind = crm_rknd ) c2(in1,idl1,ip)
  real ( kind = crm_rknd ) cc(in1,ido,ip,l1)
  real ( kind = crm_rknd ) ,intent(inout) :: ch(in2,ido,l1,ip)
  real ( kind = crm_rknd ) ch2(in2,idl1,ip)
  real ( kind = crm_rknd ) dc2
  real ( kind = crm_rknd ) dcp
  real ( kind = crm_rknd ) ds2
  real ( kind = crm_rknd ) dsp
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ic
  integer ( kind = 4 ) idij
  integer ( kind = 4 ) idp2
  integer ( kind = 4 ) ik
  integer ( kind = 4 ) im1
  integer ( kind = 4 ) im2
  integer ( kind = 4 ) ipp2
  integer ( kind = 4 ) ipph
  integer ( kind = 4 ) is
  integer ( kind = 4 ) j
  integer ( kind = 4 ) j2
  integer ( kind = 4 ) jc
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  integer ( kind = 4 ) lc
  integer ( kind = 4 ) m
  integer ( kind = 4 ) m1
  integer ( kind = 4 ) m1d
  integer ( kind = 4 ) m2
  integer ( kind = 4 ) m2s
  integer ( kind = 4 ) nbd
  real ( kind = crm_rknd ) tpi
  real ( kind = crm_rknd ) wa(ido)

  m1d = ( m - 1 ) * im1 + 1
  m2s = 1 - im2
  tpi = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 )
  arg = tpi / real ( ip, kind = crm_rknd )
  dcp = cos ( arg )
  dsp = sin ( arg )
  ipph = ( ip + 1 ) / 2
  ipp2 = ip + 2
  idp2 = ido + 2
  nbd = ( ido - 1 ) / 2

  if ( ido == 1 ) then
    do ik = 1, idl1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        c2(m1,ik,1) = ch2(m2,ik,1)
      end do
    end do

  else

    do ik = 1, idl1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch2(m2,ik,1) = c2(m1,ik,1)
      end do
    end do

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

    if ( l1 < nbd ) then

      is = -ido

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

    else

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

    end if

    if ( nbd < l1 ) then

      do j = 2, ipph
        jc = ipp2 - j
        do i = 3, ido, 2
          do k = 1, l1
            m2 = m2s
            do m1 = 1, m1d, im1
              m2 = m2 + im2
              c1(m1,i-1,k,j)  = ch(m2,i-1,k,j)  + ch(m2,i-1,k,jc)
              c1(m1,i-1,k,jc) = ch(m2,i,k,j)    - ch(m2,i,k,jc)
              c1(m1,i,k,j)    = ch(m2,i,k,j)    + ch(m2,i,k,jc)
              c1(m1,i,k,jc)   = ch(m2,i-1,k,jc) - ch(m2,i-1,k,j)
            end do
          end do
        end do
      end do

    else

      do j = 2, ipph
        jc = ipp2 - j
        do k = 1, l1
          do i = 3, ido, 2
            m2 = m2s
            do m1 = 1, m1d, im1
              m2 = m2 + im2
              c1(m1,i-1,k,j)  = ch(m2,i-1,k,j)  + ch(m2,i-1,k,jc)
              c1(m1,i-1,k,jc) = ch(m2,i,k,j)    - ch(m2,i,k,jc)
              c1(m1,i,k,j)    = ch(m2,i,k,j)    + ch(m2,i,k,jc)
              c1(m1,i,k,jc)   = ch(m2,i-1,k,jc) - ch(m2,i-1,k,j)
            end do
          end do
        end do
      end do

    end if

  end if

  do j = 2, ipph
    jc = ipp2 - j
    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        c1(m1,1,k,j)  = ch(m2,1,k,j)  + ch(m2,1,k,jc)
        c1(m1,1,k,jc) = ch(m2,1,k,jc) - ch(m2,1,k,j)
      end do
    end do
  end do

  ar1 = 1.0E+00
  ai1 = 0.0E+00

  do l = 2, ipph

    lc = ipp2 - l
    ar1h = dcp * ar1 - dsp * ai1
    ai1 =  dcp * ai1 + dsp * ar1
    ar1 = ar1h

    do ik = 1, idl1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch2(m2,ik,l)  = c2(m1,ik,1) + ar1 * c2(m1,ik,2)
        ch2(m2,ik,lc) =               ai1 * c2(m1,ik,ip)
      end do
    end do

    dc2 = ar1
    ds2 = ai1
    ar2 = ar1
    ai2 = ai1

    do j = 3, ipph
      jc = ipp2 - j
      ar2h = dc2 * ar2 - ds2 * ai2
      ai2  = dc2 * ai2 + ds2 * ar2
      ar2 = ar2h
      do ik = 1, idl1
        m2 = m2s
        do m1 = 1, m1d, im1
          m2 = m2 + im2
          ch2(m2,ik,l)  = ch2(m2,ik,l)  + ar2 * c2(m1,ik,j)
          ch2(m2,ik,lc) = ch2(m2,ik,lc) + ai2 * c2(m1,ik,jc)
        end do
      end do
    end do

  end do

  do j = 2, ipph
    do ik = 1, idl1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        ch2(m2,ik,1) = ch2(m2,ik,1) + c2(m1,ik,j)
      end do
    end do
  end do

  if ( ido < l1 ) then

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

  else

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

  end if

  do j = 2, ipph
    jc = ipp2 - j
    j2 = j + j
    do k = 1, l1
      m2 = m2s
      do m1 = 1, m1d, im1
        m2 = m2 + im2
        cc(m1,ido,j2-2,k) = ch(m2,1,k,j)
        cc(m1,1,j2-1,k) = ch(m2,1,k,jc)
      end do
    end do
  end do

  if ( ido == 1 ) then
    return
  end if

  if ( nbd < l1 ) then

    do j = 2, ipph
      jc = ipp2 - j
      j2 = j + j
      do i = 3, ido, 2
        ic = idp2 - i
        do k = 1, l1
          m2 = m2s
          do m1 = 1, m1d, im1
            m2 = m2 + im2
            cc(m1,i-1,j2-1,k)  = ch(m2,i-1,k,j) + ch(m2,i-1,k,jc)
            cc(m1,ic-1,j2-2,k) = ch(m2,i-1,k,j) - ch(m2,i-1,k,jc)
            cc(m1,i,j2-1,k)    = ch(m2,i,k,j)   + ch(m2,i,k,jc)
            cc(m1,ic,j2-2,k)   = ch(m2,i,k,jc)  - ch(m2,i,k,j)
          end do
        end do
      end do
    end do

  else

    do j = 2, ipph
      jc = ipp2 - j
      j2 = j + j
      do k = 1, l1
        do i = 3, ido, 2
          ic = idp2 - i
          m2 = m2s
          do m1 = 1, m1d, im1
            m2 = m2 + im2
            cc(m1,i-1,j2-1,k)  = ch(m2,i-1,k,j) + ch(m2,i-1,k,jc)
            cc(m1,ic-1,j2-2,k) = ch(m2,i-1,k,j) - ch(m2,i-1,k,jc)
            cc(m1,i,j2-1,k)    = ch(m2,i,k,j)   + ch(m2,i,k,jc)
            cc(m1,ic,j2-2,k)   = ch(m2,i,k,jc)  - ch(m2,i,k,j)
          end do
        end do
      end do
    end do

  end if

  return
end
!========================================================================================
! MRFTB1 is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mrftb1 (m,im,n,inc,c,ch,wa,fac)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) inc
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real ( kind = crm_rknd ) ,intent(inout) :: c(inc,*)
  real ( kind = crm_rknd ) ,intent(inout) :: ch(m,*)
  real ( kind = crm_rknd ) fac(15)
  real ( kind = crm_rknd ) half
  real ( kind = crm_rknd ) halfm
  integer ( kind = 4 ) i
  integer ( kind = 4 ) idl1
  integer ( kind = 4 ) ido
  integer ( kind = 4 ) im
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) iw
  integer ( kind = 4 ) ix2
  integer ( kind = 4 ) ix3
  integer ( kind = 4 ) ix4
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k1
  integer ( kind = 4 ) l1
  integer ( kind = 4 ) l2
  integer ( kind = 4 ) m2
  integer ( kind = 4 ) modn
  integer ( kind = 4 ) na
  integer ( kind = 4 ) nf
  integer ( kind = 4 ) nl
  real ( kind = crm_rknd ) wa(n)

  nf = int ( fac(2) )
  na = 0

  do k1 = 1, nf

    ip = int ( fac(k1+2) )
    na = 1 - na

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

  end do

  half = 0.5E+00
  halfm = -0.5E+00
  modn = mod ( n, 2 )
  nl = n - 2
  if ( modn /= 0 ) then
    nl = n - 1
  end if

  if ( na == 0 ) then

    do j = 2, nl, 2
      m2 = 1 - im
      do i = 1, m
        m2 = m2 + im
        c(m2,j) = half * c(m2,j)
        c(m2,j+1) = halfm * c(m2,j+1)
      end do
    end do

  else

    m2 = 1 - im

    do i = 1, m
      m2 = m2 + im
      ch(i,1) = c(m2,1)
      ch(i,n) = c(m2,n)
    end do

    do j = 2, nl, 2
      m2 = 1 - im
      do i = 1, m
        m2 = m2 + im
        ch(i,j) = half * c(m2,j)
        ch(i,j+1) = halfm * c(m2,j+1)
      end do
    end do

  end if

  l1 = 1
  iw = 1

  do k1 = 1, nf

    ip = int ( fac(k1+2) )
    l2 = ip * l1
    ido = n / l2
    idl1 = ido * l1

    if ( ip == 4 ) then

      ix2 = iw + ido
      ix3 = ix2 + ido

      if ( na == 0 ) then
        call mradb4 ( m, ido, l1, c, im, inc, ch, 1, m, wa(iw), wa(ix2), &
          wa(ix3) )
      else
        call mradb4 ( m, ido, l1, ch, 1, m, c, im, inc, wa(iw), wa(ix2), &
          wa(ix3) )
      end if

      na = 1 - na

    else if ( ip == 2 ) then

      if ( na == 0 ) then
        call mradb2 ( m, ido, l1, c, im, inc, ch, 1, m, wa(iw) )
      else
        call mradb2 ( m, ido, l1, ch, 1, m, c, im, inc, wa(iw) )
      end if

      na = 1 - na

    else if ( ip == 3 ) then

      ix2 = iw + ido

      if ( na == 0 ) then
        call mradb3 ( m, ido, l1, c, im, inc, ch, 1, m, wa(iw), wa(ix2) )
      else
        call mradb3 ( m, ido, l1, ch, 1, m, c, im, inc, wa(iw), wa(ix2) )
      end if

      na = 1 - na

    else if ( ip == 5 ) then

      ix2 = iw + ido
      ix3 = ix2 + ido
      ix4 = ix3 + ido

      if ( na == 0 ) then
        call mradb5 ( m, ido, l1, c, im, inc, ch, 1, m, wa(iw), wa(ix2), &
          wa(ix3), wa(ix4) )
      else
        call mradb5 ( m, ido, l1, ch, 1, m, c, im, inc, wa(iw), wa(ix2), &
          wa(ix3), wa(ix4) )
      end if

      na = 1 - na

    else

      if ( na == 0 ) then
        call mradbg ( m, ido, ip, l1, idl1, c, c, c, im, inc, ch, ch, 1, &
          m, wa(iw) )
      else
        call mradbg ( m, ido, ip, l1, idl1, ch, ch, ch, 1, m, c, c, im, &
          inc, wa(iw) )
      end if

      if ( ido == 1 ) then
        na = 1 - na
      end if

    end if

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

  end do

  return
end
!========================================================================================
! MRFTF1 is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mrftf1 (m,im,n,inc,c,ch,wa,fac)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) m
  integer ( kind = 4 ) im
  integer ( kind = 4 ) n
  integer ( kind = 4 ) inc
  real ( kind = crm_rknd ) ,intent(inout) :: c(inc,*)
  real ( kind = crm_rknd ) ,intent(inout) :: ch(m,*)
  real ( kind = crm_rknd ) ,intent(inout) :: wa(n)
  real ( kind = crm_rknd ) ,intent(inout) :: fac(15)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) idl1
  integer ( kind = 4 ) ido
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) iw
  integer ( kind = 4 ) ix2
  integer ( kind = 4 ) ix3
  integer ( kind = 4 ) ix4
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k1
  integer ( kind = 4 ) kh
  integer ( kind = 4 ) l1
  integer ( kind = 4 ) l2
  integer ( kind = 4 ) m2
  integer ( kind = 4 ) modn
  integer ( kind = 4 ) na
  integer ( kind = 4 ) nf
  integer ( kind = 4 ) nl
  real ( kind = crm_rknd ) sn
  real ( kind = crm_rknd ) tsn
  real ( kind = crm_rknd ) tsnm

  nf = int ( fac(2) )   ! number of factors
  na = 1
  l2 = n
  iw = n

! n = nx
! m = lot  = nz

! do i = 1,15
! write(*,*)'whannah - fac(',i,') : '+fac(i)
! end do

  do k1 = 1, nf

    kh = nf - k1
    ip = int ( fac(kh+3) )
! ip = fac(nf-k1+3)
! nf = fac(2)
! fac = wsave(n+1)
! ip = fac( fac(2) - k1 + 3 )

! write(*,222) k1,ip,nf,na,(kh+3),fac(kh+3)
! 222 format('whannah - k1=',i2,' ip=',i2,' nf=',i2,' na=',i2,' fac(',i2,')=',f4.1)

    l1 = l2 / ip
    ido = n / l2
    idl1 = ido * l1
    iw = iw - ( ip - 1 ) * ido
    na = 1 - na

    if ( ip == 4 ) then

      ix2 = iw + ido
      ix3 = ix2 + ido

      if ( na == 0 ) then
        call mradf4 ( m, ido, l1, c, im, inc, ch, 1,m, wa(iw), wa(ix2), &
          wa(ix3) )
      else
        call mradf4 ( m, ido, l1, ch, 1, m, c, im, inc, wa(iw), wa(ix2), &
          wa(ix3) )
      end if

    else if ( ip == 2 ) then

      if ( na == 0 ) then
        call mradf2 ( m, ido, l1, c, im, inc, ch, 1, m, wa(iw) )
      else
        call mradf2 ( m, ido, l1, ch, 1, m, c, im, inc, wa(iw) )
      end if

    else if ( ip == 3 ) then

      ix2 = iw + ido

      if ( na == 0 ) then
        call mradf3 ( m, ido, l1, c, im, inc, ch, 1, m, wa(iw), wa(ix2) )
      else
        call mradf3 ( m, ido, l1, ch, 1, m, c, im, inc, wa(iw), wa(ix2) )
      end if

    else if ( ip == 5 ) then

      ix2 = iw + ido
      ix3 = ix2 + ido
      ix4 = ix3 + ido

      if ( na == 0 ) then
        call mradf5 ( m, ido, l1, c, im, inc, ch, 1, m, wa(iw), wa(ix2), &
          wa(ix3), wa(ix4) )
      else
        call mradf5 ( m, ido, l1, ch, 1, m, c, im, inc, wa(iw), wa(ix2), &
          wa(ix3), wa(ix4) )
      end if

    else

      if ( ido == 1 ) then
        na = 1 - na
      end if

      if ( na == 0 ) then
        call mradfg ( m, ido, ip, l1, idl1, c, c, c, im, inc, ch, ch, 1, &
          m, wa(iw) )
        na = 1
      else
        call mradfg ( m, ido, ip, l1, idl1, ch, ch, ch, 1, m, c, c, im, &
          inc, wa(iw) )
        na = 0
      end if

    end if

    l2 = l1

  end do

  sn  = 1.0E+00 / real ( n, kind = crm_rknd )
  tsn = 2.0E+00 / real ( n, kind = crm_rknd )
  tsnm = -tsn
  modn = mod ( n, 2 )

  if ( modn /= 0 ) then
    nl = n - 1
  else
    nl = n - 2
  end if

  if ( na == 0 ) then

    m2 = 1-im
    do i = 1, m
      m2 = m2 + im
      c(m2,1) = sn * ch(i,1)
    end do

    do j = 2, nl, 2
      m2 = 1 - im
      do i = 1, m
        m2 = m2 + im
        c(m2,j) = tsn * ch(i,j)
        c(m2,j+1) = tsnm * ch(i,j+1)
      end do
    end do

    if ( modn == 0 ) then
      m2 = 1 - im
      do i = 1, m
        m2 = m2 + im
        c(m2,n) = sn * ch(i,n)
      end do
    end if

  else

    m2 = 1-im
    do i = 1, m
      m2 = m2 + im
      c(m2,1) = sn * c(m2,1)
    end do

    do j = 2, nl, 2
      m2 = 1 - im
      do i = 1, m
        m2 = m2 + im
        c(m2,j) = tsn * c(m2,j)
        c(m2,j+1) = tsnm * c(m2,j+1)
      end do
    end do

    if ( modn == 0 ) then

      m2 = 1 - im

      do i = 1, m
        m2 = m2 + im
        c(m2,n) = sn * c(m2,n)
      end do

    end if

  end if

  return
end
!========================================================================================
! MRFTI1 is an FFTPACK5.1 auxilliary function.
!========================================================================================
subroutine mrfti1 (n,wa,fac)

  !*****************************************************************************80
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) N, the number for which factorization and 
  !    other information is needed.
  !
  !    Output, real ( kind = crm_rknd ) WA(N), trigonometric information.
  !
  !    Output, real ( kind = crm_rknd ) FAC(15), factorization information.  FAC(1) is 
  !    N, FAC(2) is NF, the number of factors, and FAC(3:NF+2) are the factors.
  !
  implicit none

  integer ( kind = 4 )     ,intent(in   ) :: n
  real ( kind = crm_rknd ) ,intent(inout) :: wa(n)
  real ( kind = crm_rknd ) ,intent(inout) :: fac(15)

  real ( kind = 8 ) arg
  real ( kind = 8 ) argh
  real ( kind = 8 ) argld
  
  real ( kind = crm_rknd ) fi
  real ( kind = 8 ) tpi
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ib
  integer ( kind = 4 ) ido
  integer ( kind = 4 ) ii
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) ipm
  integer ( kind = 4 ) is
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k1
  integer ( kind = 4 ) l1
  integer ( kind = 4 ) l2
  integer ( kind = 4 ) ld
  integer ( kind = 4 ) nf
  integer ( kind = 4 ) nfm1
  integer ( kind = 4 ) nl
  integer ( kind = 4 ) nq
  integer ( kind = 4 ) nr
  integer ( kind = 4 ) ntry
  integer ( kind = 4 ) ntryh(4)
  
  

  nl = n
  nf = 0
  j = 0

  do while ( 1 < nl )

    j = j + 1

    if ( j == 1 ) then
      ntry = 4
    else if ( j == 2 ) then
      ntry = 2
    else if ( j == 3 ) then
      ntry = 3
    else if ( j == 4 ) then
      ntry = 5
    else
      ntry = ntry + 2
    end if

    do

      nq = nl / ntry
      nr = nl - ntry * nq

      if ( nr /= 0 ) then
        exit
      end if

      nf = nf + 1
      fac(nf+2) = real ( ntry, kind = crm_rknd )
      nl = nq
!
!  If 2 is a factor, make sure it appears first in the list of factors.
!
      if ( ntry == 2 ) then
        if ( nf /= 1 ) then
          do i = 2, nf
            ib = nf - i + 2
            fac(ib+2) = fac(ib+1)
          end do
          fac(3) = 2.0E+00
        end if
      end if

    end do

  end do

  fac(1) = real ( n , kind = crm_rknd )
  fac(2) = real ( nf, kind = crm_rknd )
  tpi = 8.0D+00 * atan ( 1.0D+00 )
  argh = tpi / real ( n, kind = crm_rknd )
  is = 0
  nfm1 = nf - 1
  l1 = 1

  do k1 = 1, nfm1
    ip = int ( fac(k1+2) )
    ld = 0
    l2 = l1 * ip
    ido = n / l2
    ipm = ip - 1
    do j = 1, ipm
      ld = ld + l1
      i = is
      argld = real ( ld, kind = crm_rknd ) * argh
      fi = 0.0E+00
      do ii = 3, ido, 2
        i = i + 2
        fi = fi + 1.0E+00
        arg = fi * argld
        wa(i-1) = cos ( arg )
        wa(i) = sin ( arg )
      end do
      is = is + ido
    end do
    l1 = l2
  end do

  return
end
!========================================================================================
! XERCON checks INC, JUMP, N and LOT for consistency.
!========================================================================================
function xercon ( inc, jump, n, lot )

  !*****************************************************************************80
  !
  !  Discussion:
  !
  !    Positive integers INC, JUMP, N and LOT are "consistent" if,
  !    for any values I1 and I2 < N, and J1 and J2 < LOT,
  !
  !      I1 * INC + J1 * JUMP = I2 * INC + J2 * JUMP
  !
  !    can only occur if I1 = I2 and J1 = J2.
  !
  !    For multiple FFT's to execute correctly, INC, JUMP, N and LOT must
  !    be consistent, or else at least one array element will be
  !    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:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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 ( kind = 4 ) INC, JUMP, N, LOT, the parameters to check.
  !
  !    Output, logical XERCON, is TRUE if the parameters are consistent.
  !
  implicit none

  integer ( kind = 4 ) i
  integer ( kind = 4 ) inc
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jnew
  integer ( kind = 4 ) jump
  integer ( kind = 4 ) lcm
  integer ( kind = 4 ) lot
  integer ( kind = 4 ) n
  logical xercon

  i = inc
  j = jump

  do while ( j /= 0 )
    jnew = mod ( i, j )
    i = j
    j = jnew
  end do
!
!  LCM = least common multiple of INC and JUMP.
!
  lcm = ( inc * jump ) / i

  if ( lcm <= ( n - 1 ) * inc .and. lcm <= ( lot - 1 ) * jump ) then
    xercon = .false.
  else
    xercon = .true.
  end if

  return
end
!========================================================================================
! XERFFT is an error handler for the FFTPACK routines.
!========================================================================================
subroutine xerfft ( srname, info )

  !*****************************************************************************80
  !
  !  Discussion:
  !
  !    XERFFT is an error handler for FFTPACK version 5.1 routines.
  !    It is called by an FFTPACK 5.1 routine if an input parameter has an
  !    invalid value.  A message is printed and execution stops.
  !
  !    Installers may consider modifying the stop statement in order to
  !    call system-specific exception-handling facilities.
  !
  !  License:
  !
  !    Licensed under the GNU General Public License (GPL).
  !    Copyright (C) 1995-2004, Scientific Computing Division,
  !    University Corporation for Atmospheric Research
  !
  !  Modified:
  !
  !    31 July 2011
  !
  !  Author:
  !
  !    Original FORTRAN77 version by Paul Swarztrauber, Richard Valent.
  !    FORTRAN90 version by John Burkardt.
  !
  !  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, character ( len = * ) SRNAME, the name of the calling routine.
  !
  !    Input, integer ( kind = 4 ) INFO, an error code.  When a single invalid 
  !    parameter in the parameter list of the calling routine has been detected, 
  !    INFO is the position of that parameter.  In the case when an illegal 
  !    combination of LOT, JUMP, N, and INC has been detected, the calling 
  !    subprogram calls XERFFT with INFO = -1.
  !
  implicit none

  integer ( kind = 4 ) info
  character ( len = * ) srname

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'XERFFT - Fatal error!'

  if ( 1 <= info ) then
    write ( *, '(a,a,a,i3,a)') '  On entry to ', trim ( srname ), &
      ' parameter number ', info, ' had an illegal value.'
  else if ( info == -1 ) then
    write( *, '(a,a,a,a)') '  On entry to ', trim ( srname ), &
      ' parameters LOT, JUMP, N and INC are inconsistent.'
  else if ( info == -2 ) then
    write( *, '(a,a,a,a)') '  On entry to ', trim ( srname ), &
      ' parameter L is greater than LDIM.'
  else if ( info == -3 ) then
    write( *, '(a,a,a,a)') '  On entry to ', trim ( srname ), &
      ' parameter M is greater than MDIM.'
  else if ( info == -5 ) then
    write( *, '(a,a,a,a)') '  Within ', trim ( srname ), &
      ' input error returned by lower level routine.'
  else if ( info == -6 ) then
    write( *, '(a,a,a,a)') '  On entry to ', trim ( srname ), &
      ' parameter LDIM is less than 2*(L/2+1).'
  end if

  stop
end
!========================================================================================
!========================================================================================
end module fftpack5
