fft.f90

!MNH_LIC Copyright 1994-2023 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
! Modifications:
!  P. Wautelet 22/02/2019: replace Hollerith edit descriptor (deleted from Fortran 95 standard)
!  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
!  P. Wautelet 22/11/2022: modernization of FFts (fixed-format F77 -> free-format F95)
!-----------------------------------------------------------------
MODULE MODE_FFT

USE MODE_MSG

IMPLICIT NONE

PRIVATE

PUBLIC :: FFT991, SET99
PUBLIC :: NVECLEN

INTEGER, PARAMETER :: NVECLEN=1020*16

REAL, PARAMETER :: XSIN36 = 0.587785252292473
REAL, PARAMETER :: XSIN60 = 0.866025403784437
REAL, PARAMETER :: XSIN72 = 0.95105651629515
REAL, PARAMETER :: XQRT5  = 0.559016994374947

CONTAINS

SUBROUTINE SET99( TRIGS, IFAX, N )
  ! SUBROUTINE 'SET99' - COMPUTES FACTORS OF N & TRIGONOMETRIC
  ! FUNCTIONS REQUIRED BY FFT991

  IMPLICIT NONE

  REAL,    DIMENSION(N),  INTENT(INOUT) :: TRIGS
  INTEGER, DIMENSION(19), INTENT(INOUT) :: IFAX
  INTEGER,                INTENT(IN)    :: N

  INTEGER, PARAMETER, DIMENSION(7) :: NLFAX = [ 6, 8, 5, 4, 3, 2, 1 ]

  INTEGER :: I, IFAC, IL, K, NFAX, NU
  INTEGER, DIMENSION(10) :: JFAX
  REAL :: ANGLE, DEL

  DEL=4.0*ASIN(1.0)/REAL(N)
  DO K = 0, (N/2)-1
    ANGLE=REAL(K)*DEL
    TRIGS(2*K+1)=COS(ANGLE)
    TRIGS(2*K+2)=SIN(ANGLE)
  END DO

  ! FIND FACTORS OF N (8,6,5,4,3,2; ONLY ONE 8 ALLOWED)
  ! LOOK FOR SIXES FIRST, STORE FACTORS IN DESCENDING ORDER
  NU=N
  IFAC=6
  K=0
  IL=1

  DO
    IF ( MOD(NU,IFAC) == 0 ) THEN
      K=K+1
      JFAX(K)=IFAC
      IF ( IFAC == 8 .AND. K /= 1 ) THEN
        JFAX(1)=8
        JFAX(K)=6
      END IF
      NU=NU/IFAC
      IF ( NU == 1 ) EXIT
      IF ( IFAC /= 8 ) CYCLE
    END IF

    IL=IL+1
    IFAC=NLFAX(IL)
    IF ( IFAC <= 1 ) CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'SET99', 'illegal factor' )
  END DO

  ! NOW REVERSE ORDER OF FACTORS
  NFAX=K
  IFAX(1)=NFAX
  DO I=1,NFAX
    IFAX(NFAX+2-I)=JFAX(I)
  END DO
  IFAX(10)=N

END SUBROUTINE SET99


SUBROUTINE FFT991( A, WORK, TRIGS, IFAX, JUMP, N, ILOT, ISIGN, KSZA, KSZW, KSZT )
  USE MODE_MPPDB

  IMPLICIT NONE

      REAL,    DIMENSION(KSZA), INTENT(INOUT) :: A
      REAL,    DIMENSION(KSZW), INTENT(OUT) :: WORK
      REAL,    DIMENSION(KSZT), INTENT(IN) :: TRIGS
      INTEGER, DIMENSION(19),   INTENT(IN) :: IFAX
      INTEGER,                  INTENT(IN) :: JUMP, N, ILOT, ISIGN
      INTEGER,                  INTENT(IN) :: KSZA, KSZW, KSZT
!
!     SUBROUTINE 'FFT991' - MULTIPLE FAST REAL PERIODIC TRANSFORM
!     SUPERSEDES PREVIOUS ROUTINE 'FFT991'
!
!     REAL TRANSFORM OF LENGTH N PERFORMED BY REMOVING REDUNDANT
!     OPERATIONS FROM COMPLEX TRANSFORM OF LENGTH N
!
!     A IS THE ARRAY CONTAINING INPUT & OUTPUT DATA
!     WORK IS AN AREA OF SIZE (N+1)*MIN(ILOT,LOT)
!     TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
!     IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N
!     JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
!     N IS THE LENGTH OF THE DATA VECTORS
!     ILOT IS THE NUMBER OF DATA VECTORS
!     ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
!           = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
!
!     ORDERING OF COEFFICIENTS:
!         A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2)
!         WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED
!
!     ORDERING OF DATA:
!         X(0),X(1),X(2),...,X(N-1), 0 , 0 ; (N+2) LOCATIONS REQUIRED
!
!     VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS
!     IN PARALLEL
!
!     N MUST BE COMPOSED OF FACTORS 2,3 & 5 BUT DOES NOT HAVE TO BE EVEN
!
!     DEFINITION OF TRANSFORMS:
!     -------------------------
!
!     ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))
!         WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K)
!
!     ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N))
!               B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N))
!
  INTEGER :: NFAX, NX, NBLOX, NVEX
  INTEGER :: I, J, IA, ISTART, ILA, IGO
  INTEGER :: NB, K, IFAC, IERR
  INTEGER :: IBASE, JBASE
  INTEGER :: II, JJ, IX, IZ
  INTEGER :: I0

  IF (MPPDB_INITIALIZED) THEN
    !Check all IN arrays
    CALL MPPDB_CHECK( TRIGS, "FFT991 beg:TRIGS" )
    CALL MPPDB_CHECK( IFAX, "FFT991 beg:IFAX" )
    !Check all INOUT arrays
    CALL MPPDB_CHECK( A, "FFT991 beg:A" )
  END IF

!$acc data present( A, WORK )

  ! Initialisation of WORK useful to compare results with MPPDB_CHECK (otherwise all values are not set by FFT991
  WORK(:) = 0.

#if 0
  !PW: Original version: but in that case, intent of TRIGS and IFAX are not correct
  IF(IFAX(10).NE.N) CALL SET99(TRIGS,IFAX,N)
#else
  IF(IFAX(10).NE.N) CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'FFT991', 'TRIGS and IFAX not initialised (set99 not yet called)' )
#endif
  NFAX=IFAX(1)
  NX=N+1
  IF (MOD(N,2).EQ.1) NX=N
  NBLOX=1+(ILOT-1)/NVECLEN
  NVEX=ILOT-(NBLOX-1)*NVECLEN

  IF ( ISIGN == 1 ) THEN
!
!   ISIGN=+1, SPECTRAL TO GRIDPOINT TRANSFORM
!   -----------------------------------------
    ISTART=1
    DO NB=1,NBLOX
      IA=ISTART
!$acc kernels
!CDIR NODEP
!*vocl loop,novrec
!$acc loop independent private( I )
      DO J=1,NVEX
        I = ISTART + ( J - 1 ) * JUMP
        A(I+1)=0.5*A(I)
      END DO
!$acc end kernels
      IF ( MOD(N,2) == 0 ) THEN
!$acc kernels
        I0 = ISTART + N
!CDIR NODEP
!*vocl loop,novrec
!$acc loop independent private( I )
        DO J=1,NVEX
          I = I0 + ( J - 1 ) * JUMP
          A(I)=0.5*A(I)
        END DO
!$acc end kernels
      END IF
      IA=ISTART+1
      ILA=1
      IGO=+1
!
      DO K=1,NFAX
        IFAC=IFAX(K+1)
        IERR=-1
        IF ( IGO == 1 ) THEN
          CALL RPASSM(A(IA:),A(IA+ILA:),WORK(1:),WORK(IFAC*ILA+1:), &
              TRIGS(:),                                             &
              JUMP,NX,NVEX,N,IFAC,ILA,IERR,                         &
              SIZE(A(IA:)),SIZE(A(IA+ILA:)),SIZE(WORK(1:)),         &
              SIZE(WORK(IFAC*ILA+1:)),SIZE(TRIGS(:)))
        ELSE
          CALL RPASSM(WORK(1:),WORK(ILA+1:),A(IA:),A(IA+IFAC*ILA:), &
             TRIGS(:),                                              &
             NX,JUMP,NVEX,N,IFAC,ILA,IERR,                          &
             SIZE(WORK(1:)),SIZE(WORK(ILA+1:)),SIZE(A(IA:)),        &
             SIZE(A(IA+IFAC*ILA:)),SIZE(TRIGS(:)))
        END IF
        IF (IERR.NE.0) GO TO 500
        ILA=IFAC*ILA
        IGO=-IGO
        IA=ISTART
      END DO
!
!     IF NECESSARY, COPY RESULTS BACK TO A
!     ------------------------------------
      IF ( MOD(NFAX,2) == 1 ) THEN
!$acc kernels
!CDIR NODEP
!*vocl loop,novrec
!$acc loop independent private( IBASE, JBASE )
        DO JJ = 1, NVEX
          IBASE = 1  + (JJ-1) * NX
          JBASE = IA + (JJ-1) * JUMP
          A(JBASE:JBASE+N-1) = WORK(IBASE:IBASE+N-1)
        END DO
!$acc end kernels
      END IF
!
!     FILL IN ZEROS AT END
!     --------------------
!$acc kernels
      A(ISTART+N::JUMP) = 0.0
      A(ISTART+N+1::JUMP) = 0.0
!$acc end kernels
!
      ISTART=ISTART+NVEX*JUMP
      NVEX=NVECLEN
    END DO

  ELSE IF ( ISIGN == -1 ) THEN
!
!     ISIGN=-1, GRIDPOINT TO SPECTRAL TRANSFORM
!     -----------------------------------------
    ISTART=1
    DO NB=1,NBLOX
      IA=ISTART
      ILA=N
      IGO=+1
!
      DO K=1,NFAX
        IFAC=IFAX(NFAX+2-K)
        ILA=ILA/IFAC
        IERR=-1
        IF ( IGO == 1 ) THEN
          CALL QPASSM(A(IA:),A(IA+IFAC*ILA:),WORK(1:),WORK(ILA+1:), &
              TRIGS(:),                                             &
              JUMP,NX,NVEX,N,IFAC,ILA,IERR,                         &
              SIZE(A(IA:)),SIZE(A(IA+IFAC*ILA:)),SIZE(WORK(1:)),    &
              SIZE(WORK(ILA+1:)),SIZE(TRIGS(:)))
        ELSE
          CALL QPASSM(WORK(1:),WORK(IFAC*ILA+1:),A(IA:),A(IA+ILA:), &
              TRIGS(:),                                             &
              NX,JUMP,NVEX,N,IFAC,ILA,IERR,                         &
              SIZE(WORK(1:)),SIZE(WORK(IFAC*ILA+1:)),SIZE(A(IA:)),  &
              SIZE(A(IA+ILA:)),SIZE(TRIGS(:)))
        END IF
        IF (IERR.NE.0) GO TO 500
        IGO=-IGO
        IA=ISTART+1
      END DO
!
!     IF NECESSARY, COPY RESULTS BACK TO A
!     ------------------------------------
      IF ( MOD(NFAX,2) == 1 ) THEN
!$acc kernels
        !CDIR NODEP
!*vocl loop,novrec
!$acc loop independent private( IBASE, JBASE )
        DO JJ = 1, NVEX
          IBASE = 1  + (JJ-1) * NX
          JBASE = IA + (JJ-1) * JUMP
          A(JBASE:JBASE+N-1) = WORK(IBASE:IBASE+N-1)
        END DO
!$acc end kernels
      END IF
!
!     SHIFT A(0) & FILL IN ZERO IMAG PARTS
!     ------------------------------------
!$acc kernels
!CDIR NODEP
!*vocl loop,novrec
!$acc loop independent private( I )
      DO J=1,NVEX
        IX  = ISTART + ( J - 1 ) * JUMP
        A(IX)=A(IX+1)
        A(IX+1)=0.0
      END DO
!$acc end kernels
      IF ( MOD(N,2) == 0 ) THEN
!$acc kernels
        A(ISTART+N+1::JUMP) = 0.0
!$acc end kernels
      END IF
!
      ISTART=ISTART+NVEX*JUMP
      NVEX=NVECLEN
    END DO
  END IF
!
!     ERROR MESSAGES
!     --------------
  500 CONTINUE
      IF ( IERR /= 0 ) THEN
      IF(IERR.EQ.1) THEN
  520 FORMAT('VECTOR LENGTH =',I4,', GREATER THAN NVECLEN')
      ELSEIF(IERR.EQ.2) THEN
      WRITE(6,540) IFAC
  540 FORMAT( 'FACTOR =',I3,', NOT CATERED FOR')
      ELSE
      WRITE(6,560) IFAC
  560 FORMAT('FACTOR =',I3,', ONLY CATERED FOR IF ILA*IFAC=N')
      ENDIF
      END IF

!$acc end data

  IF (MPPDB_INITIALIZED) THEN
    !Check all INOUT arrays
    CALL MPPDB_CHECK( A, "FFT991 end:A" )
    !Check all OUT arrays
    CALL MPPDB_CHECK( WORK, "FFT991 end:WORK" )
  END IF
END SUBROUTINE FFT991



SUBROUTINE RPASSM(A,B,C,D,TRIGS,INC3,INC4,ILOT,N,IFAC,ILA,IERR,KSZ1,KSZ2,KSZ3,KSZ4,KSZ5)

  IMPLICIT NONE

  REAL, DIMENSION(KSZ1), INTENT(IN)    :: A
  REAL, DIMENSION(KSZ2), INTENT(IN)    :: B
  REAL, DIMENSION(KSZ3), INTENT(INOUT) :: C
  REAL, DIMENSION(KSZ4), INTENT(INOUT) :: D
  REAL, DIMENSION(KSZ5), INTENT(IN)    :: TRIGS
  INTEGER,               INTENT(IN)    :: INC3, INC4, ILOT, N, IFAC, ILA
  INTEGER,               INTENT(OUT)   :: IERR
  INTEGER,               INTENT(IN)    :: KSZ1,KSZ2,KSZ3,KSZ4,KSZ5
!
!     SUBROUTINE 'RPASSM' - PERFORMS ONE PASS THROUGH DATA AS PART
!     OF MULTIPLE REAL FFT (FOURIER SYNTHESIS) ROUTINE
!
!     A IS FIRST REAL INPUT VECTOR
!         EQUIVALENCE B(1) WITH A (ILA+1)
!     C IS FIRST REAL OUTPUT VECTOR
!         EQUIVALENCE D(1) WITH C(IFAC*ILA+1)
!     TRIGS IS A PRECALCULATED LIST OF SINES & COSINES
!     INC3 IS THE INCREMENT BETWEEN INPUT VECTORS A
!     INC4 IS THE INCREMENT BETWEEN OUTPUT VECTORS C
!     ILOT IS THE NUMBER OF VECTORS
!     N IS THE LENGTH OF THE VECTORS
!     IFAC IS THE CURRENT FACTOR OF N
!     ILA IS THE PRODUCT OF PREVIOUS FACTORS
!     IERR IS AN ERROR INDICATOR:
!              0 - PASS COMPLETED WITHOUT ERROR
!              1 - ILOT GREATER THAN NVECLEN
!              2 - IFAC NOT CATERED FOR
!              3 - IFAC ONLY CATERED FOR IF ILA=N/IFAC
!
!-----------------------------------------------------------------------
!
  INTEGER :: M, IINK, JINK, JUMP, KSTOP
  INTEGER :: IBAD, IBASE, JBASE, IGO
  INTEGER :: I, J, IA, IB, JA, JB
  INTEGER :: IL, IJK, K, KB, IC, JC, KC
  INTEGER :: ID, JD, KD, IE, JE, KE, IF, JF, KF, JG, JH
  INTEGER :: IA0, IB0, IC0, ID0, IE0, IF0, JBASE0
  REAL :: A10, A11, A20, A21, B10, B11, B20, B21
  REAL :: C1, C2, C3, C4, C5, S1, S2, S3, S4, S5
  REAL :: QQRT5, SIN45, SSIN36, SSIN45, SSIN60, SSIN72


!$acc data present( A, B, C, D ) copyin( TRIGS )

!acc kernels
      M=N/IFAC
      IINK=ILA
      JINK=ILA
      JUMP=(IFAC-1)*JINK
      KSTOP=(N-IFAC)/(2*IFAC)

!acc end kernels
      IBAD=1
      IF (ILOT.GT.NVECLEN) GO TO 910
      IBASE=0
      JBASE=0
      IGO=IFAC-1
      IF (IGO.EQ.7) IGO=6
      IBAD=2
      IF (IGO.LT.1.OR.IGO.GT.6) GO TO 910

      SELECT CASE ( IFAC )

      CASE ( 2 )
!
!     CODING FOR FACTOR 2
!     -------------------
      IA=1
      IB=IA+2*M-ILA
      JA=1
      JB=JA+JINK

      IF ( ILA /= M ) THEN

!$acc kernels
!$acc loop independent
        DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
          DO IJK=1,ILOT
            I = IBASE + IL - 1 + (IJK - 1 ) * INC3
            J = JBASE + IL - 1 + (IJK - 1 ) * INC4
            C(JA+J)=A(IA+I)+A(IB+I)
            C(JB+J)=A(IA+I)-A(IB+I)
          END DO
        END DO
        IBASE=IBASE+ILA
        JBASE=JBASE+ILA
!$acc end kernels
        IA=IA+IINK
        IINK=2*IINK
        IB=IB-IINK
        IBASE=0
        JBASE=JBASE+JUMP
        JUMP=2*JUMP+JINK
        IF ( IA /= IB ) THEN
!$acc kernels
          JBASE0 = JBASE
          IA0 = IA
          IB0 = IB
!$acc loop independent private( KB, C1, S1, JBASE, IA, IB )
          DO K=ILA,KSTOP,ILA
            KB=K+K
            C1=TRIGS(KB+1)
            S1=TRIGS(KB+2)
            JBASE = JBASE0 + (K-ILA)/ILA * (ILA + JUMP )
            IA = IA0 + (K-ILA)/ILA * IINK
            IB = IB0 - (K-ILA)/ILA * IINK
!$acc loop independent
            DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
              DO IJK=1,ILOT
                I =         IL - 1 + (IJK - 1 ) * INC3
                J = JBASE + IL - 1 + (IJK - 1 ) * INC4
                C(JA+J)=A(IA+I)+A(IB+I)
                D(JA+J)=B(IA+I)-B(IB+I)
                C(JB+J)=C1*(A(IA+I)-A(IB+I))-S1*(B(IA+I)+B(IB+I))
                D(JB+J)=S1*(A(IA+I)-A(IB+I))+C1*(B(IA+I)+B(IB+I))
              END DO
            END DO
          END DO
          JBASE = JBASE0 + ((KSTOP-ILA)/ILA+1) * ( ILA + JUMP )
          IA = IA0 + ((KSTOP-ILA)/ILA+1) * IINK
          IB = IB0 - ((KSTOP-ILA)/ILA+1) * IINK
!$acc end kernels
        END IF
        IF ( IA <= IB ) THEN
!$acc kernels
          IBASE=0
!$acc loop independent
          DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
            DO IJK=1,ILOT
              I = IBASE + IL - 1 + (IJK - 1 ) * INC3
              J = JBASE + IL - 1 + (IJK - 1 ) * INC4
              C(JA+J)=A(IA+I)
              C(JB+J)=-B(IA+I)
            END DO
          END DO
          IBASE=IBASE+ILA
          JBASE=JBASE+ILA
!$acc end kernels
        END IF

      ELSE

!$acc kernels
!$acc loop independent
        DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
          DO IJK=1,ILOT
            I = IBASE + IL - 1 + (IJK - 1 ) * INC3
            J = JBASE + IL - 1 + (IJK - 1 ) * INC4
            C(JA+J)=2.0*(A(IA+I)+A(IB+I))
            C(JB+J)=2.0*(A(IA+I)-A(IB+I))
          END DO
        END DO
        IBASE=IBASE+ILA
        JBASE=JBASE+ILA
!$acc end kernels
      END IF
      IBAD=0


      CASE ( 3 )
!
!     CODING FOR FACTOR 3
!     -------------------
      IA=1
      IB=IA+2*M-ILA
      IC=IB
      JA=1
      JB=JA+JINK
      JC=JB+JINK

      IF ( ILA /= M ) THEN

!$acc kernels
!$acc loop independent
        DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
          DO IJK=1,ILOT
            I = IBASE + IL - 1 + (IJK - 1 ) * INC3
            J = JBASE + IL - 1 + (IJK - 1 ) * INC4
            C(JA+J)=A(IA+I)+A(IB+I)
            C(JB+J)=(A(IA+I)-0.5*A(IB+I))-(XSIN60*(B(IB+I)))
            C(JC+J)=(A(IA+I)-0.5*A(IB+I))+(XSIN60*(B(IB+I)))
          END DO
        END DO
        IBASE=IBASE+ILA
        JBASE=JBASE+ILA
!$acc end kernels
        IA=IA+IINK
        IINK=2*IINK
        IB=IB+IINK
        IC=IC-IINK
        JBASE=JBASE+JUMP
        JUMP=2*JUMP+JINK
        IF ( IA /= IC ) THEN
!$acc kernels
          JBASE0 = JBASE
          IA0 = IA
          IB0 = IB
          IC0 = IC
!$acc loop independent private( KB, KC, C1, S1, C2, S2, C3, S3, JBASE, IA, IB, IC )
          DO K=ILA,KSTOP,ILA
            KB=K+K
            KC=KB+KB
            C1=TRIGS(KB+1)
            S1=TRIGS(KB+2)
            C2=TRIGS(KC+1)
            S2=TRIGS(KC+2)
            JBASE = JBASE0 + (K-ILA)/ILA * (ILA + JUMP )
            IA = IA0 + (K-ILA)/ILA * IINK
            IB = IB0 + (K-ILA)/ILA * IINK
            IC = IC0 - (K-ILA)/ILA * IINK
!$acc loop independent
            DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
              DO IJK=1,ILOT
                I =         IL - 1 + (IJK - 1 ) * INC3
                J = JBASE + IL - 1 + (IJK - 1 ) * INC4
                C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
                D(JA+J)=B(IA+I)+(B(IB+I)-B(IC+I))
                C(JB+J)=                                                      &
                    C1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(XSIN60*(B(IB+I)+B(IC+ &
                  I))))                                                       &
                  -S1*((B(IA+I)-0.5*(B(IB+I)-B(IC+I)))+(XSIN60*(A(IB+I)-A(IC+ &
                  I))))
                D(JB+J)=                                                      &
                    S1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(XSIN60*(B(IB+I)+B(IC+ &
                  I))))                                                       &
                  +C1*((B(IA+I)-0.5*(B(IB+I)-B(IC+I)))+(XSIN60*(A(IB+I)-A(IC+ &
                  I))))
                C(JC+J)=                                                      &
                    C2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(XSIN60*(B(IB+I)+B(IC+ &
                  I))))                                                       &
                  -S2*((B(IA+I)-0.5*(B(IB+I)-B(IC+I)))-(XSIN60*(A(IB+I)-A(IC+ &
                  I))))
                D(JC+J)=                                                      &
                    S2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(XSIN60*(B(IB+I)+B(IC+ &
                  I))))                                                       &
                  +C2*((B(IA+I)-0.5*(B(IB+I)-B(IC+I)))-(XSIN60*(A(IB+I)-A(IC+ &
                  I))))
              END DO
            END DO
          END DO
          JBASE = JBASE0 + ((KSTOP-ILA)/ILA+1) * ( ILA + JUMP )
          IA = IA0 + ((KSTOP-ILA)/ILA+1) * IINK
          IB = IB0 + ((KSTOP-ILA)/ILA+1) * IINK
          IC = IC0 - ((KSTOP-ILA)/ILA+1) * IINK
!$acc end kernels
        END IF

        IF ( IA <= IC ) THEN
!$acc kernels
          IBASE=0
!$acc loop independent
          DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
            DO IJK=1,ILOT
              I = IBASE + IL - 1 + (IJK - 1 ) * INC3
              J = JBASE + IL - 1 + (IJK - 1 ) * INC4
              C(JA+J)=A(IA+I)+A(IB+I)
              C(JB+J)=(0.5*A(IA+I)-A(IB+I))-(XSIN60*B(IA+I))
              C(JC+J)=-(0.5*A(IA+I)-A(IB+I))-(XSIN60*B(IA+I))
            END DO
          END DO
          IBASE = IBASE + ILA
          JBASE = JBASE + ILA
!$acc end kernels
        END IF

      ELSE

!$acc kernels
        SSIN60=2.0*XSIN60
!$acc loop independent
        DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
          DO IJK=1,ILOT
            I = IBASE + IL - 1 + (IJK - 1 ) * INC3
            J = JBASE + IL - 1 + (IJK - 1 ) * INC4
            C(JA+J)=2.0*(A(IA+I)+A(IB+I))
            C(JB+J)=(2.0*A(IA+I)-A(IB+I))-(SSIN60*B(IB+I))
            C(JC+J)=(2.0*A(IA+I)-A(IB+I))+(SSIN60*B(IB+I))
          END DO
        END DO
        IBASE = IBASE + ILA
        JBASE = JBASE + ILA
!$acc end kernels
      END IF

      IBAD=0


      CASE ( 4 )
!
!     CODING FOR FACTOR 4
!     -------------------
      IA=1
      IB=IA+2*M-ILA
      IC=IB+2*M
      ID=IB
      JA=1
      JB=JA+JINK
      JC=JB+JINK
      JD=JC+JINK

      IF ( ILA /=  M) THEN

!$acc kernels
!$acc loop independent
        DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
          DO IJK=1,ILOT
            I = IBASE + IL - 1 + (IJK - 1 ) * INC3
            J = JBASE + IL - 1 + (IJK - 1 ) * INC4
            C(JA+J)=(A(IA+I)+A(IC+I))+A(IB+I)
            C(JB+J)=(A(IA+I)-A(IC+I))-B(IB+I)
            C(JC+J)=(A(IA+I)+A(IC+I))-A(IB+I)
            C(JD+J)=(A(IA+I)-A(IC+I))+B(IB+I)
          END DO
        END DO
        IBASE = IBASE + ILA
        JBASE = JBASE + ILA
!$acc end kernels
        IA=IA+IINK
        IINK=2*IINK
        IB=IB+IINK
        IC=IC-IINK
        ID=ID-IINK
        JBASE=JBASE+JUMP
        JUMP=2*JUMP+JINK
        IF ( IB /=  IC ) THEN
!$acc kernels
          JBASE0 = JBASE
          IA0 = IA
          IB0 = IB
          IC0 = IC
          ID0 = ID
!$acc loop independent private( KB, KC, KD, C1, S1, C2, S2, C3, S3, JBASE, IA, IB, IC, ID )
          DO K=ILA,KSTOP,ILA
            KB=K+K
            KC=KB+KB
            KD=KC+KB
            C1=TRIGS(KB+1)
            S1=TRIGS(KB+2)
            C2=TRIGS(KC+1)
            S2=TRIGS(KC+2)
            C3=TRIGS(KD+1)
            S3=TRIGS(KD+2)
            JBASE = JBASE0 + (K-ILA)/ILA * (ILA + JUMP )
            IA = IA0 + (K-ILA)/ILA * IINK
            IB = IB0 + (K-ILA)/ILA * IINK
            IC = IC0 - (K-ILA)/ILA * IINK
            ID = ID0 - (K-ILA)/ILA * IINK
!$acc loop independent
            DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
              DO IJK=1,ILOT
                I =         IL - 1 + (IJK - 1 ) * INC3
                J = JBASE + IL - 1 + (IJK - 1 ) * INC4
                C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))
                D(JA+J)=(B(IA+I)-B(IC+I))+(B(IB+I)-B(ID+I))
                C(JC+J)=                                     &
                    C2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) &
                   -S2*((B(IA+I)-B(IC+I))-(B(IB+I)-B(ID+I)))
                D(JC+J)=                                     &
                    S2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) &
                   +C2*((B(IA+I)-B(IC+I))-(B(IB+I)-B(ID+I)))
                C(JB+J)=                                     &
                    C1*((A(IA+I)-A(IC+I))-(B(IB+I)+B(ID+I))) &
                   -S1*((B(IA+I)+B(IC+I))+(A(IB+I)-A(ID+I)))
                D(JB+J)=                                     &
                    S1*((A(IA+I)-A(IC+I))-(B(IB+I)+B(ID+I))) &
                   +C1*((B(IA+I)+B(IC+I))+(A(IB+I)-A(ID+I)))
                C(JD+J)=                                     &
                    C3*((A(IA+I)-A(IC+I))+(B(IB+I)+B(ID+I))) &
                   -S3*((B(IA+I)+B(IC+I))-(A(IB+I)-A(ID+I)))
                D(JD+J)=                                     &
                    S3*((A(IA+I)-A(IC+I))+(B(IB+I)+B(ID+I))) &
                   +C3*((B(IA+I)+B(IC+I))-(A(IB+I)-A(ID+I)))
              END DO
            END DO
          END DO
          JBASE = JBASE0 + ((KSTOP-ILA)/ILA+1) * ( ILA + JUMP )
          IA = IA0 + ((KSTOP-ILA)/ILA+1) * IINK
          IB = IB0 + ((KSTOP-ILA)/ILA+1) * IINK
          IC = IC0 - ((KSTOP-ILA)/ILA+1) * IINK
          ID = ID0 - ((KSTOP-ILA)/ILA+1) * IINK
!$acc end kernels
        END IF

        IF ( IB <= IC ) THEN
!$acc kernels
          IBASE=0
          SIN45=SQRT(0.5)
!$acc loop independent
          DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
            DO IJK=1,ILOT
              I = IBASE + IL - 1 + (IJK - 1 ) * INC3
              J = JBASE + IL - 1 + (IJK - 1 ) * INC4
              C(JA+J)=A(IA+I)+A(IB+I)
              C(JB+J)=SIN45*((A(IA+I)-A(IB+I))-(B(IA+I)+B(IB+I)))
              C(JC+J)=B(IB+I)-B(IA+I)
              C(JD+J)=-SIN45*((A(IA+I)-A(IB+I))+(B(IA+I)+B(IB+I)))
            END DO
          END DO
          IBASE = IBASE + ILA
          JBASE = JBASE + ILA
!$acc end kernels
        END IF

      ELSE ! ILA == M

!$acc kernels
!$acc loop independent
        DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
          DO IJK=1,ILOT
            I = IBASE + IL - 1 + (IJK - 1 ) * INC3
            J = JBASE + IL - 1 + (IJK - 1 ) * INC4
            C(JA+J)=2.0*((A(IA+I)+A(IC+I))+A(IB+I))
            C(JB+J)=2.0*((A(IA+I)-A(IC+I))-B(IB+I))
            C(JC+J)=2.0*((A(IA+I)+A(IC+I))-A(IB+I))
            C(JD+J)=2.0*((A(IA+I)-A(IC+I))+B(IB+I))
          END DO
        END DO
        IBASE = IBASE + ILA
        JBASE = JBASE + ILA
!$acc end kernels
      END IF

      IBAD=0


      CASE ( 5 )
!
!     CODING FOR FACTOR 5
!     -------------------
      IA=1
      IB=IA+2*M-ILA
      IC=IB+2*M
      ID=IC
      IE=IB
      JA=1
      JB=JA+JINK
      JC=JB+JINK
      JD=JC+JINK
      JE=JD+JINK

      IF ( ILA /=  M ) THEN

!$acc kernels
!$acc loop independent
        DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
          DO IJK=1,ILOT
            I = IBASE + IL - 1 + (IJK - 1 ) * INC3
            J = JBASE + IL - 1 + (IJK - 1 ) * INC4
            C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
            C(JB+J)=((A(IA+I)-0.25*(A(IB+I)+A(IC+I)))+XQRT5*(A(IB+I)-A(IC+ &
              I)))                                                        &
                -(XSIN72*B(IB+I)+XSIN36*B(IC+I))
            C(JC+J)=((A(IA+I)-0.25*(A(IB+I)+A(IC+I)))-XQRT5*(A(IB+I)-A(IC+ &
              I)))                                                        &
                -(XSIN36*B(IB+I)-XSIN72*B(IC+I))
            C(JD+J)=((A(IA+I)-0.25*(A(IB+I)+A(IC+I)))-XQRT5*(A(IB+I)-A(IC+ &
              I)))                                                        &
                +(XSIN36*B(IB+I)-XSIN72*B(IC+I))
            C(JE+J)=((A(IA+I)-0.25*(A(IB+I)+A(IC+I)))+XQRT5*(A(IB+I)-A(IC+ &
              I)))                                                        &
                +(XSIN72*B(IB+I)+XSIN36*B(IC+I))
          END DO
        END DO
        IBASE = IBASE + ILA
        JBASE = JBASE + ILA
!$acc end kernels
        IA=IA+IINK
        IINK=2*IINK
        IB=IB+IINK
        IC=IC+IINK
        ID=ID-IINK
        IE=IE-IINK
        JBASE=JBASE+JUMP
        JUMP=2*JUMP+JINK
        IF ( IB /= ID ) THEN
!$acc kernels
          JBASE0 = JBASE
          IA0 = IA
          IB0 = IB
          IC0 = IC
          ID0 = ID
          IE0 = IE
!$acc loop independent private( KB, KC, KD, KE, C1, S1, C2, S2, C3, S3, C4, S4, JBASE, IA, IB, IC, ID, IE )
          DO K=ILA,KSTOP,ILA
            KB=K+K
            KC=KB+KB
            KD=KC+KB
            KE=KD+KB
            C1=TRIGS(KB+1)
            S1=TRIGS(KB+2)
            C2=TRIGS(KC+1)
            S2=TRIGS(KC+2)
            C3=TRIGS(KD+1)
            S3=TRIGS(KD+2)
            C4=TRIGS(KE+1)
            S4=TRIGS(KE+2)
            JBASE = JBASE0 + (K-ILA)/ILA * (ILA + JUMP )
            IA = IA0 + (K-ILA)/ILA * IINK
            IB = IB0 + (K-ILA)/ILA * IINK
            IC = IC0 + (K-ILA)/ILA * IINK
            ID = ID0 - (K-ILA)/ILA * IINK
            IE = IE0 - (K-ILA)/ILA * IINK
!$acc loop independent
            DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J, A10, A11, A20, A21, B10, B11, B20, B21 )
              DO IJK=1,ILOT
                I =         IL - 1 + (IJK - 1 ) * INC3
                J = JBASE + IL - 1 + (IJK - 1 ) * INC4

                A10=(A(IA+I)-0.25*((A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I)))) &
                    +XQRT5*((A(IB+I)+A(IE+I))-(A(IC+I)+A(ID+I)))
                A20=(A(IA+I)-0.25*((A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I)))) &
                    -XQRT5*((A(IB+I)+A(IE+I))-(A(IC+I)+A(ID+I)))
                B10=(B(IA+I)-0.25*((B(IB+I)-B(IE+I))+(B(IC+I)-B(ID+I)))) &
                    +XQRT5*((B(IB+I)-B(IE+I))-(B(IC+I)-B(ID+I)))
                B20=(B(IA+I)-0.25*((B(IB+I)-B(IE+I))+(B(IC+I)-B(ID+I)))) &
                    -XQRT5*((B(IB+I)-B(IE+I))-(B(IC+I)-B(ID+I)))
                A11=XSIN72*(B(IB+I)+B(IE+I))+XSIN36*(B(IC+I)+B(ID+I))
                A21=XSIN36*(B(IB+I)+B(IE+I))-XSIN72*(B(IC+I)+B(ID+I))
                B11=XSIN72*(A(IB+I)-A(IE+I))+XSIN36*(A(IC+I)-A(ID+I))
                B21=XSIN36*(A(IB+I)-A(IE+I))-XSIN72*(A(IC+I)-A(ID+I))

                C(JA+J)=A(IA+I)+((A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I)))
                D(JA+J)=B(IA+I)+((B(IB+I)-B(IE+I))+(B(IC+I)-B(ID+I)))
                C(JB+J)=C1*(A10-A11)-S1*(B10+B11)
                D(JB+J)=S1*(A10-A11)+C1*(B10+B11)
                C(JE+J)=C4*(A10+A11)-S4*(B10-B11)
                D(JE+J)=S4*(A10+A11)+C4*(B10-B11)
                C(JC+J)=C2*(A20-A21)-S2*(B20+B21)
                D(JC+J)=S2*(A20-A21)+C2*(B20+B21)
                C(JD+J)=C3*(A20+A21)-S3*(B20-B21)
                D(JD+J)=S3*(A20+A21)+C3*(B20-B21)
              END DO
            END DO
          END DO
          JBASE = JBASE0 + ((KSTOP-ILA)/ILA+1) * ( ILA + JUMP )
          IA = IA0 + ((KSTOP-ILA)/ILA+1) * IINK
          IB = IB0 + ((KSTOP-ILA)/ILA+1) * IINK
          IC = IC0 + ((KSTOP-ILA)/ILA+1) * IINK
          ID = ID0 - ((KSTOP-ILA)/ILA+1) * IINK
          IE = IE0 - ((KSTOP-ILA)/ILA+1) * IINK
!$acc end kernels
        END IF

        IF ( IB <= ID ) THEN
!$acc kernels
          IBASE=0
!$acc loop independent
          DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
            DO IJK=1,ILOT
              I = IBASE + IL - 1 + (IJK - 1 ) * INC3
              J = JBASE + IL - 1 + (IJK - 1 ) * INC4
              C(JA+J)=(A(IA+I)+A(IB+I))+A(IC+I)
              C(JB+J)=(XQRT5*(A(IA+I)-A(IB+I))+(0.25*(A(IA+I)+A(IB+I))-A(IC+ &
                I)))                                                        &
                  -(XSIN36*B(IA+I)+XSIN72*B(IB+I))
              C(JE+J)=-(XQRT5*(A(IA+I)-A(IB+I))+(0.25*(A(IA+I)+A(IB+I))-A(IC+ &
                I)))                                                         &
                  -(XSIN36*B(IA+I)+XSIN72*B(IB+I))
              C(JC+J)=(XQRT5*(A(IA+I)-A(IB+I))-(0.25*(A(IA+I)+A(IB+I))-A(IC+ &
                I)))                                                        &
                  -(XSIN72*B(IA+I)-XSIN36*B(IB+I))
              C(JD+J)=-(XQRT5*(A(IA+I)-A(IB+I))-(0.25*(A(IA+I)+A(IB+I))-A(IC+ &
                I)))                                                         &
                  -(XSIN72*B(IA+I)-XSIN36*B(IB+I))
            END DO
          END DO
          IBASE = IBASE + ILA
          JBASE = JBASE + ILA
!$acc end kernels
        END IF

      ELSE

!$acc kernels
        QQRT5=2.0*XQRT5
        SSIN36=2.0*XSIN36
        SSIN72=2.0*XSIN72
!$acc loop independent
        DO IL=1,ILA
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( I, J )
          DO IJK=1,ILOT
            I = IBASE + IL - 1 + (IJK - 1 ) * INC3
            J = JBASE + IL - 1 + (IJK - 1 ) * INC4
            C(JA+J)=2.0*(A(IA+I)+(A(IB+I)+A(IC+I)))
            C(JB+J)=(2.0*(A(IA+I)-0.25*(A(IB+I)+A(IC+I)))                 &
                +QQRT5*(A(IB+I)-A(IC+I)))-(SSIN72*B(IB+I)+SSIN36*B(IC+I))