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( PTRIGS, KFAX, KN )
  ! SUBROUTINE 'SET99' - COMPUTES FACTORS OF KN & TRIGONOMETRIC
  ! FUNCTIONS REQUIRED BY FFT991

  IMPLICIT NONE

  REAL,    DIMENSION(KN), INTENT(INOUT) :: PTRIGS
  INTEGER, DIMENSION(19), INTENT(INOUT) :: KFAX
  INTEGER,                INTENT(IN)    :: KN

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

  INTEGER :: II, IFAC, IL, IK, IFAX, IU
  INTEGER, DIMENSION(10) :: IJFAX
  REAL :: ZANGLE, ZDEL

  ZDEL=4.0*ASIN(1.0)/REAL(KN)
  DO IK = 0, (KN/2)-1
    ZANGLE=REAL(IK)*ZDEL
    PTRIGS(2*IK+1)=COS(ZANGLE)
    PTRIGS(2*IK+2)=SIN(ZANGLE)
  END DO

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

  DO
    IF ( MOD(IU,IFAC) == 0 ) THEN
      IK=IK+1
      IJFAX(IK)=IFAC
      IF ( IFAC == 8 .AND. IK /= 1 ) THEN
        IJFAX(1)=8
        IJFAX(IK)=6
      END IF
      IU=IU/IFAC
      IF ( IU == 1 ) EXIT
      IF ( IFAC /= 8 ) CYCLE
    END IF

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

  ! NOW REVERSE ORDER OF FACTORS
  IFAX=IK
  KFAX(1)=IFAX
  DO II=1,IFAX
    KFAX(IFAX+2-II)=IJFAX(II)
  END DO
  KFAX(10)=KN

END SUBROUTINE SET99


SUBROUTINE FFT991( PA, PWORK, PTRIGS, PFAX, KJUMP, KN, KLOT, KSIGN )
  USE MODE_MPPDB

  IMPLICIT NONE

  REAL,    DIMENSION(:),  INTENT(INOUT) :: PA
  REAL,    DIMENSION(:,:,:), CONTIGUOUS, TARGET, INTENT(OUT) :: PWORK
  REAL,    DIMENSION(:),  INTENT(IN)    :: PTRIGS
  INTEGER, DIMENSION(19), INTENT(IN)    :: PFAX
  INTEGER,                INTENT(IN)    :: KJUMP, KN, KLOT, KSIGN
!
!     SUBROUTINE 'FFT991' - MULTIPLE FAST REAL PERIODIC TRANSFORM
!     SUPERSEDES PREVIOUS ROUTINE 'FFT991'
!
!     REAL TRANSFORM OF LENGTH KN PERFORMED BY REMOVING REDUNDANT
!     OPERATIONS FROM COMPLEX TRANSFORM OF LENGTH KN
!
!     PA IS THE ARRAY CONTAINING INPUT & OUTPUT DATA
!     PWORK IS AN AREA OF SIZE (KN+1)*MIN(KLOT,LOT)
!     PTRIGS IS PA PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
!     KFAX IS PA PREVIOUSLY PREPARED LIST OF FACTORS OF KN
!     JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
!     KN IS THE LENGTH OF THE DATA VECTORS
!     KLOT IS THE NUMBER OF DATA VECTORS
!     KSIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
!           = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
!
!     ORDERING OF COEFFICIENTS:
!         PA(0),PB(0),PA(1),PB(1),PA(2),PB(2),...,PA(KN/2),PB(KN/2)
!         WHERE PB(0)=PB(KN/2)=0; (KN+2) LOCATIONS REQUIRED
!
!     ORDERING OF DATA:
!         X(0),X(1),X(2),...,X(KN-1), 0 , 0 ; (KN+2) LOCATIONS REQUIRED
!
!     VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS
!     IN PARALLEL
!
!     KN MUST BE COMPOSED OF FACTORS 2,3 & 5 BUT DOES NOT HAVE TO BE EVEN
!
!     DEFINITION OF TRANSFORMS:
!     -------------------------
!
!     KSIGN=+1: X(J)=SUM(K=0,...,KN-1)(PC(K)*EXP(2*I*J*K*PI/KN))
!         WHERE PC(K)=PA(K)+I*PB(K) AND PC(KN-K)=PA(K)-I*PB(K)
!
!     KSIGN=-1: PA(K)=(1/KN)*SUM(J=0,...,KN-1)(X(J)*COS(2*J*K*PI/KN))
!               PB(K)=-(1/KN)*SUM(J=0,...,KN-1)(X(J)*SIN(2*J*K*PI/KN))
!
  INTEGER :: IFAX, INX, INBLOX, INVEX
  INTEGER :: I, J, IA, ISTART, ILA, IGO
  INTEGER :: INB, IK, IFAC, IERR
  INTEGER :: IIBASE, IJBASE
  INTEGER :: II, JJ, IX, IZ
  INTEGER :: I0
  REAL, DIMENSION(:), POINTER :: ZWORK

  IF (MPPDB_INITIALIZED) THEN
    !Check all IN arrays
    CALL MPPDB_CHECK( PTRIGS, "FFT991 beg:PTRIGS" )
    CALL MPPDB_CHECK( PFAX, "FFT991 beg:KFAX" )
    !Check all INOUT arrays
    CALL MPPDB_CHECK( PA, "FFT991 beg:PA" )
  END IF

ZWORK(1 : SIZE(PWORK)) => PWORK(:,:,:)

!$acc data present( PA, ZWORK )

  ! Initialisation of ZWORK useful to compare results with MPPDB_CHECK (otherwise all values are not set by FFT991
  ZWORK(:) = 0.
!$acc update device( ZWORK )

#if 0
  !PW: Original version: but in that case, intent of PTRIGS and KFAX are not correct
  IF(PFAX(10).NE.KN) CALL SET99(PTRIGS,KFAX,KN)
#else
  IF(PFAX(10).NE.KN) CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'FFT991', 'PTRIGS and KFAX not initialised (set99 not yet called)' )
#endif
  IFAX=PFAX(1)
  INX=KN+1
  IF (MOD(KN,2).EQ.1) INX=KN
  INBLOX=1+(KLOT-1)/NVECLEN
  INVEX=KLOT-(INBLOX-1)*NVECLEN

  IF ( KSIGN == 1 ) THEN
!
!   KSIGN=+1, SPECTRAL TO GRIDPOINT TRANSFORM
!   -----------------------------------------
    ISTART=1
    DO INB=1,INBLOX
      IA=ISTART
!$acc kernels
!CDIR NODEP
!*vocl loop,novrec
!$acc loop independent private( I )
      DO J=1,INVEX
        I = ISTART + ( J - 1 ) * KJUMP
        PA(I+1)=0.5*PA(I)
      END DO
!$acc end kernels
      IF ( MOD(KN,2) == 0 ) THEN
!$acc kernels
        I0 = ISTART + KN
!CDIR NODEP
!*vocl loop,novrec
!$acc loop independent private( I )
        DO J=1,INVEX
          I = I0 + ( J - 1 ) * KJUMP
          PA(I)=0.5*PA(I)
        END DO
!$acc end kernels
      END IF
      IA=ISTART+1
      ILA=1
      IGO=+1
!
      DO IK=1,IFAX
        IFAC=PFAX(IK+1)
        IERR=-1
        IF ( IGO == 1 ) THEN
          CALL RPASSM( PA(IA:), PA(IA+ILA:), ZWORK(:), ZWORK(IFAC*ILA+1:), PTRIGS(:), &
                       KJUMP, INX, INVEX, KN, IFAC, ILA, IERR )
        ELSE
          CALL RPASSM( ZWORK(:), ZWORK(ILA+1:), PA(IA:), PA(IA+IFAC*ILA:), PTRIGS(:), &
                       INX, KJUMP, INVEX, KN, IFAC, ILA, IERR )
        END IF
        IF (IERR.NE.0) GO TO 500
        ILA=IFAC*ILA
        IGO=-IGO
        IA=ISTART
      END DO
!
!     IF NECESSARY, COPY RESULTS BACK TO PA
!     -------------------------------------
      IF ( MOD(IFAX,2) == 1 ) THEN
!$acc kernels
!CDIR NODEP
!*vocl loop,novrec
!$acc loop independent private( IIBASE, IJBASE )
        DO JJ = 1, INVEX
          IIBASE = 1  + (JJ-1) * INX
          IJBASE = IA + (JJ-1) * KJUMP
          PA(IJBASE:IJBASE+KN-1) = ZWORK(IIBASE:IIBASE+KN-1)
        END DO
!$acc end kernels
      END IF
!
!     FILL IN ZEROS AT END
!     --------------------
!$acc kernels
      PA(ISTART+KN::KJUMP) = 0.0
      PA(ISTART+KN+1::KJUMP) = 0.0
!$acc end kernels
!
      ISTART=ISTART+INVEX*KJUMP
      INVEX=NVECLEN
    END DO

  ELSE IF ( KSIGN == -1 ) THEN
!
!     KSIGN=-1, GRIDPOINT TO SPECTRAL TRANSFORM
!     -----------------------------------------
    ISTART=1
    DO INB=1,INBLOX
      IA=ISTART
      ILA=KN
      IGO=+1
!
      DO IK=1,IFAX
        IFAC=PFAX(IFAX+2-IK)
        ILA=ILA/IFAC
        IERR=-1
        IF ( IGO == 1 ) THEN
          CALL QPASSM( PA(IA:), PA(IA+IFAC*ILA:), ZWORK(:), ZWORK(ILA+1:), PTRIGS(:), &
                       KJUMP, INX, INVEX, KN, IFAC, ILA, IERR )
        ELSE
          CALL QPASSM( ZWORK(:), ZWORK(IFAC*ILA+1:), PA(IA:), PA(IA+ILA:), PTRIGS(:), &
                       INX, KJUMP, INVEX, KN, IFAC, ILA, IERR )
        END IF
        IF (IERR.NE.0) GO TO 500
        IGO=-IGO
        IA=ISTART+1
      END DO
!
!     IF NECESSARY, COPY RESULTS BACK TO PA
!     -------------------------------------
      IF ( MOD(IFAX,2) == 1 ) THEN
!$acc kernels
        !CDIR NODEP
!*vocl loop,novrec
!$acc loop independent private( IIBASE, IJBASE )
        DO JJ = 1, INVEX
          IIBASE = 1  + (JJ-1) * INX
          IJBASE = IA + (JJ-1) * KJUMP
          PA(IJBASE:IJBASE+KN-1) = ZWORK(IIBASE:IIBASE+KN-1)
        END DO
!$acc end kernels
      END IF
!
!     SHIFT PA(0) & FILL IN ZERO IMAG PARTS
!     -------------------------------------
!$acc kernels
!CDIR NODEP
!*vocl loop,novrec
!$acc loop independent private( I )
      DO J=1,INVEX
        IX  = ISTART + ( J - 1 ) * KJUMP
        PA(IX)=PA(IX+1)
        PA(IX+1)=0.0
      END DO
!$acc end kernels
      IF ( MOD(KN,2) == 0 ) THEN
!$acc kernels
        PA(ISTART+KN+1::KJUMP) = 0.0
!$acc end kernels
      END IF
!
      ISTART=ISTART+INVEX*KJUMP
      INVEX=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=KN')
      ENDIF
      END IF

!$acc end data

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



SUBROUTINE RPASSM( PA, PB, PC, PD, PTRIGS, KINC3, KINC4, KLOT, KN, KFAC, KLA, KERR )

  IMPLICIT NONE

  REAL, DIMENSION(:), INTENT(IN)    :: PA
  REAL, DIMENSION(:), INTENT(IN)    :: PB
  REAL, DIMENSION(:), INTENT(INOUT) :: PC
  REAL, DIMENSION(:), INTENT(INOUT) :: PD
  REAL, DIMENSION(:), INTENT(IN)    :: PTRIGS
  INTEGER,            INTENT(IN)    :: KINC3, KINC4, KLOT, KN, KFAC, KLA
  INTEGER,            INTENT(OUT)   :: KERR
!
!     SUBROUTINE 'RPASSM' - PERFORMS ONE PASS THROUGH DATA AS PART
!     OF MULTIPLE REAL FFT (FOURIER SYNTHESIS) ROUTINE
!
!     KA IS FIRST REAL INPUT VECTOR
!         EQUIVALENCE PB(1) WITH PA (KLA+1)
!     PC IS FIRST REAL OUTPUT VECTOR
!         EQUIVALENCE PD(1) WITH PC(KFAC*KLA+1)
!     PTRIGS IS PA PRECALCULATED LIST OF SINES & COSINES
!     KINC3 IS THE INCREMENT BETWEEN INPUT VECTORS PA
!     KINC4 IS THE INCREMENT BETWEEN OUTPUT VECTORS PC
!     KLOT IS THE NUMBER OF VECTORS
!     KN IS THE LENGTH OF THE VECTORS
!     KFAC IS THE CURRENT FACTOR OF KN
!     KLA IS THE PRODUCT OF PREVIOUS FACTORS
!     KERR IS AN ERROR INDICATOR:
!              0 - PASS COMPLETED WITHOUT ERROR
!              1 - KLOT GREATER THAN NVECLEN
!              2 - KFAC NOT CATERED FOR
!              3 - KFAC ONLY CATERED FOR IF KLA=KN/KFAC
!
!-----------------------------------------------------------------------
!
  INTEGER :: IM, IINK, JINK, JUMP, ISTOP
  INTEGER :: IBAD, IIBASE, IJBASE, IGO
  INTEGER :: II, IJ, IIA, IIB, IJA, IJB
  INTEGER :: IIL, IJK, IIK, IKB, IIC, IJC, IKC
  INTEGER :: IID, IJD, IKD, IIE, IJE, IKE, IIF, IJF, IKF, IJG, IJH
  INTEGER :: IIA0, IIB0, IIC0, IID0, IIE0, IIF0, IJBASE0
  REAL :: ZA10, ZA11, ZA20, ZA21, ZB10, ZB11, ZB20, ZB21
  REAL :: ZC1, ZC2, ZC3, ZC4, ZC5, ZS1, ZS2, ZS3, ZS4, ZS5
  REAL :: ZQQRT5, ZSIN45, ZSSIN36, ZSSIN45, ZSSIN60, ZSSIN72


!$acc data present( PA, PB, PC, PD ) copyin( PTRIGS )

!acc kernels
      IM=KN/KFAC
      IINK=KLA
      JINK=KLA
      JUMP=(KFAC-1)*JINK
      ISTOP=(KN-KFAC)/(2*KFAC)

!acc end kernels
      IBAD=1
      IF (KLOT.GT.NVECLEN) GO TO 910
      IIBASE=0
      IJBASE=0
      IGO=KFAC-1
      IF (IGO.EQ.7) IGO=6
      IBAD=2
      IF (IGO.LT.1.OR.IGO.GT.6) GO TO 910

      SELECT CASE ( KFAC )

      CASE ( 2 )
!
!     CODING FOR FACTOR 2
!     -------------------
      IIA=1
      IIB=IIA+2*IM-KLA
      IJA=1
      IJB=IJA+JINK

      IF ( KLA /= IM ) THEN

!$acc kernels
!$acc loop independent
        DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
          DO IIL=0,KLA-1
            II = IIBASE + IJK * KINC3 + IIL
            IJ = IJBASE + IJK * KINC4 + IIL
            PC(IJA+IJ)=PA(IIA+II)+PA(IIB+II)
            PC(IJB+IJ)=PA(IIA+II)-PA(IIB+II)
          END DO
        END DO
        IIBASE=IIBASE+KLA
        IJBASE=IJBASE+KLA
!$acc end kernels
        IIA=IIA+IINK
        IINK=2*IINK
        IIB=IIB-IINK
        IIBASE=0
        IJBASE=IJBASE+JUMP
        JUMP=2*JUMP+JINK
        IF ( IIA /= IIB ) THEN
!$acc kernels
          IJBASE0 = IJBASE
          IIA0 = IIA
          IIB0 = IIB
!$acc loop independent private( IKB, ZC1, ZS1, IJBASE, IIA, IIB )
          DO IIK=KLA,ISTOP,KLA
            IKB=IIK+IIK
            ZC1=PTRIGS(IKB+1)
            ZS1=PTRIGS(IKB+2)
            IJBASE = IJBASE0 + (IIK-KLA)/KLA * (KLA + JUMP )
            IIA = IIA0 + (IIK-KLA)/KLA * IINK
            IIB = IIB0 - (IIK-KLA)/KLA * IINK
!$acc loop independent
            DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
              DO IIL=0,KLA-1
                II =          IJK * KINC3 + IIL
                IJ = IJBASE + IJK * KINC4 + IIL
                PC(IJA+IJ)=PA(IIA+II)+PA(IIB+II)
                PD(IJA+IJ)=PB(IIA+II)-PB(IIB+II)
                PC(IJB+IJ)=ZC1*(PA(IIA+II)-PA(IIB+II))-ZS1*(PB(IIA+II)+PB(IIB+II))
                PD(IJB+IJ)=ZS1*(PA(IIA+II)-PA(IIB+II))+ZC1*(PB(IIA+II)+PB(IIB+II))
              END DO
            END DO
          END DO
          IJBASE = IJBASE0 + ((ISTOP-KLA)/KLA+1) * ( KLA + JUMP )
          IIA = IIA0 + ((ISTOP-KLA)/KLA+1) * IINK
          IIB = IIB0 - ((ISTOP-KLA)/KLA+1) * IINK
!$acc end kernels
        END IF
        IF ( IIA <= IIB ) THEN
!$acc kernels
          IIBASE=0
!$acc loop independent
          DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
            DO IIL=0,KLA-1
              II = IIBASE + IJK * KINC3 + IIL
              IJ = IJBASE + IJK * KINC4 + IIL
              PC(IJA+IJ)=PA(IIA+II)
              PC(IJB+IJ)=-PB(IIA+II)
            END DO
          END DO
          IIBASE=IIBASE+KLA
          IJBASE=IJBASE+KLA
!$acc end kernels
        END IF

      ELSE

!$acc kernels
!$acc loop independent
        DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
          DO IIL=0,KLA-1
            II = IIBASE + IJK * KINC3 + IIL
            IJ = IJBASE + IJK * KINC4 + IIL
            PC(IJA+IJ)=2.0*(PA(IIA+II)+PA(IIB+II))
            PC(IJB+IJ)=2.0*(PA(IIA+II)-PA(IIB+II))
          END DO
        END DO
        IIBASE=IIBASE+KLA
        IJBASE=IJBASE+KLA
!$acc end kernels
      END IF
      IBAD=0


      CASE ( 3 )
!
!     CODING FOR FACTOR 3
!     -------------------
      IIA=1
      IIB=IIA+2*IM-KLA
      IIC=IIB
      IJA=1
      IJB=IJA+JINK
      IJC=IJB+JINK

      IF ( KLA /= IM ) THEN

!$acc kernels
!$acc loop independent
        DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
          DO IIL=0,KLA-1
            II = IIBASE + IJK * KINC3 + IIL
            IJ = IJBASE + IJK * KINC4 + IIL
            PC(IJA+IJ)=PA(IIA+II)+PA(IIB+II)
            PC(IJB+IJ)=(PA(IIA+II)-0.5*PA(IIB+II))-(XSIN60*(PB(IIB+II)))
            PC(IJC+IJ)=(PA(IIA+II)-0.5*PA(IIB+II))+(XSIN60*(PB(IIB+II)))
          END DO
        END DO
        IIBASE=IIBASE+KLA
        IJBASE=IJBASE+KLA
!$acc end kernels
        IIA=IIA+IINK
        IINK=2*IINK
        IIB=IIB+IINK
        IIC=IIC-IINK
        IJBASE=IJBASE+JUMP
        JUMP=2*JUMP+JINK
        IF ( IIA /= IIC ) THEN
!$acc kernels
          IJBASE0 = IJBASE
          IIA0 = IIA
          IIB0 = IIB
          IIC0 = IIC
!$acc loop independent private( IKB, IKC, ZC1, ZS1, ZC2, ZS2, ZC3, ZS3, IJBASE, IIA, IIB, IIC )
          DO IIK=KLA,ISTOP,KLA
            IKB=IIK+IIK
            IKC=IKB+IKB
            ZC1=PTRIGS(IKB+1)
            ZS1=PTRIGS(IKB+2)
            ZC2=PTRIGS(IKC+1)
            ZS2=PTRIGS(IKC+2)
            IJBASE = IJBASE0 + (IIK-KLA)/KLA * (KLA + JUMP )
            IIA = IIA0 + (IIK-KLA)/KLA * IINK
            IIB = IIB0 + (IIK-KLA)/KLA * IINK
            IIC = IIC0 - (IIK-KLA)/KLA * IINK
!$acc loop independent
            DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
              DO IIL=0,KLA-1
                II =          IJK * KINC3 + IIL
                IJ = IJBASE + IJK * KINC4 + IIL
                PC(IJA+IJ)=PA(IIA+II)+(PA(IIB+II)+PA(IIC+II))
                PD(IJA+IJ)=PB(IIA+II)+(PB(IIB+II)-PB(IIC+II))
                PC(IJB+IJ)=                                                      &
                    ZC1*((PA(IIA+II)-0.5*(PA(IIB+II)+PA(IIC+II)))-(XSIN60*(PB(IIB+II)+PB(IIC+ &
                  II))))                                                       &
                  -ZS1*((PB(IIA+II)-0.5*(PB(IIB+II)-PB(IIC+II)))+(XSIN60*(PA(IIB+II)-PA(IIC+ &
                  II))))
                PD(IJB+IJ)=                                                      &
                    ZS1*((PA(IIA+II)-0.5*(PA(IIB+II)+PA(IIC+II)))-(XSIN60*(PB(IIB+II)+PB(IIC+ &
                  II))))                                                       &
                  +ZC1*((PB(IIA+II)-0.5*(PB(IIB+II)-PB(IIC+II)))+(XSIN60*(PA(IIB+II)-PA(IIC+ &
                  II))))
                PC(IJC+IJ)=                                                      &
                    ZC2*((PA(IIA+II)-0.5*(PA(IIB+II)+PA(IIC+II)))+(XSIN60*(PB(IIB+II)+PB(IIC+ &
                  II))))                                                       &
                  -ZS2*((PB(IIA+II)-0.5*(PB(IIB+II)-PB(IIC+II)))-(XSIN60*(PA(IIB+II)-PA(IIC+ &
                  II))))
                PD(IJC+IJ)=                                                      &
                    ZS2*((PA(IIA+II)-0.5*(PA(IIB+II)+PA(IIC+II)))+(XSIN60*(PB(IIB+II)+PB(IIC+ &
                  II))))                                                       &
                  +ZC2*((PB(IIA+II)-0.5*(PB(IIB+II)-PB(IIC+II)))-(XSIN60*(PA(IIB+II)-PA(IIC+ &
                  II))))
              END DO
            END DO
          END DO
          IJBASE = IJBASE0 + ((ISTOP-KLA)/KLA+1) * ( KLA + JUMP )
          IIA = IIA0 + ((ISTOP-KLA)/KLA+1) * IINK
          IIB = IIB0 + ((ISTOP-KLA)/KLA+1) * IINK
          IIC = IIC0 - ((ISTOP-KLA)/KLA+1) * IINK
!$acc end kernels
        END IF

        IF ( IIA <= IIC ) THEN
!$acc kernels
          IIBASE=0
!$acc loop independent
          DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
            DO IIL=0,KLA-1
              II = IIBASE + IJK * KINC3 + IIL
              IJ = IJBASE + IJK * KINC4 + IIL
              PC(IJA+IJ)=PA(IIA+II)+PA(IIB+II)
              PC(IJB+IJ)=(0.5*PA(IIA+II)-PA(IIB+II))-(XSIN60*PB(IIA+II))
              PC(IJC+IJ)=-(0.5*PA(IIA+II)-PA(IIB+II))-(XSIN60*PB(IIA+II))
            END DO
          END DO
          IIBASE = IIBASE + KLA
          IJBASE = IJBASE + KLA
!$acc end kernels
        END IF

      ELSE

!$acc kernels
        ZSSIN60=2.0*XSIN60
!$acc loop independent
        DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
          DO IIL=0,KLA-1
            II = IIBASE + IJK * KINC3 + IIL
            IJ = IJBASE + IJK * KINC4 + IIL
            PC(IJA+IJ)=2.0*(PA(IIA+II)+PA(IIB+II))
            PC(IJB+IJ)=(2.0*PA(IIA+II)-PA(IIB+II))-(ZSSIN60*PB(IIB+II))
            PC(IJC+IJ)=(2.0*PA(IIA+II)-PA(IIB+II))+(ZSSIN60*PB(IIB+II))
          END DO
        END DO
        IIBASE = IIBASE + KLA
        IJBASE = IJBASE + KLA
!$acc end kernels
      END IF

      IBAD=0


      CASE ( 4 )
!
!     CODING FOR FACTOR 4
!     -------------------
      IIA=1
      IIB=IIA+2*IM-KLA
      IIC=IIB+2*IM
      IID=IIB
      IJA=1
      IJB=IJA+JINK
      IJC=IJB+JINK
      IJD=IJC+JINK

      IF ( KLA /=  IM) THEN

!$acc kernels
!$acc loop independent
        DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
          DO IIL=0,KLA-1
            II = IIBASE + IJK * KINC3 + IIL
            IJ = IJBASE + IJK * KINC4 + IIL
            PC(IJA+IJ)=(PA(IIA+II)+PA(IIC+II))+PA(IIB+II)
            PC(IJB+IJ)=(PA(IIA+II)-PA(IIC+II))-PB(IIB+II)
            PC(IJC+IJ)=(PA(IIA+II)+PA(IIC+II))-PA(IIB+II)
            PC(IJD+IJ)=(PA(IIA+II)-PA(IIC+II))+PB(IIB+II)
          END DO
        END DO
        IIBASE = IIBASE + KLA
        IJBASE = IJBASE + KLA
!$acc end kernels
        IIA=IIA+IINK
        IINK=2*IINK
        IIB=IIB+IINK
        IIC=IIC-IINK
        IID=IID-IINK
        IJBASE=IJBASE+JUMP
        JUMP=2*JUMP+JINK
        IF ( IIB /=  IIC ) THEN
!$acc kernels
          IJBASE0 = IJBASE
          IIA0 = IIA
          IIB0 = IIB
          IIC0 = IIC
          IID0 = IID
!$acc loop independent private( IKB, IKC, IKD, ZC1, ZS1, ZC2, ZS2, ZC3, ZS3, IJBASE, IIA, IIB, IIC, IID )
          DO IIK=KLA,ISTOP,KLA
            IKB=IIK+IIK
            IKC=IKB+IKB
            IKD=IKC+IKB
            ZC1=PTRIGS(IKB+1)
            ZS1=PTRIGS(IKB+2)
            ZC2=PTRIGS(IKC+1)
            ZS2=PTRIGS(IKC+2)
            ZC3=PTRIGS(IKD+1)
            ZS3=PTRIGS(IKD+2)
            IJBASE = IJBASE0 + (IIK-KLA)/KLA * (KLA + JUMP )
            IIA = IIA0 + (IIK-KLA)/KLA * IINK
            IIB = IIB0 + (IIK-KLA)/KLA * IINK
            IIC = IIC0 - (IIK-KLA)/KLA * IINK
            IID = IID0 - (IIK-KLA)/KLA * IINK
!$acc loop independent
            DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
              DO IIL=0,KLA-1
                II =          IJK * KINC3 + IIL
                IJ = IJBASE + IJK * KINC4 + IIL
                PC(IJA+IJ)=(PA(IIA+II)+PA(IIC+II))+(PA(IIB+II)+PA(IID+II))
                PD(IJA+IJ)=(PB(IIA+II)-PB(IIC+II))+(PB(IIB+II)-PB(IID+II))
                PC(IJC+IJ)=                                     &
                    ZC2*((PA(IIA+II)+PA(IIC+II))-(PA(IIB+II)+PA(IID+II))) &
                   -ZS2*((PB(IIA+II)-PB(IIC+II))-(PB(IIB+II)-PB(IID+II)))
                PD(IJC+IJ)=                                     &
                    ZS2*((PA(IIA+II)+PA(IIC+II))-(PA(IIB+II)+PA(IID+II))) &
                   +ZC2*((PB(IIA+II)-PB(IIC+II))-(PB(IIB+II)-PB(IID+II)))
                PC(IJB+IJ)=                                     &
                    ZC1*((PA(IIA+II)-PA(IIC+II))-(PB(IIB+II)+PB(IID+II))) &
                   -ZS1*((PB(IIA+II)+PB(IIC+II))+(PA(IIB+II)-PA(IID+II)))
                PD(IJB+IJ)=                                     &
                    ZS1*((PA(IIA+II)-PA(IIC+II))-(PB(IIB+II)+PB(IID+II))) &
                   +ZC1*((PB(IIA+II)+PB(IIC+II))+(PA(IIB+II)-PA(IID+II)))
                PC(IJD+IJ)=                                     &
                    ZC3*((PA(IIA+II)-PA(IIC+II))+(PB(IIB+II)+PB(IID+II))) &
                   -ZS3*((PB(IIA+II)+PB(IIC+II))-(PA(IIB+II)-PA(IID+II)))
                PD(IJD+IJ)=                                     &
                    ZS3*((PA(IIA+II)-PA(IIC+II))+(PB(IIB+II)+PB(IID+II))) &
                   +ZC3*((PB(IIA+II)+PB(IIC+II))-(PA(IIB+II)-PA(IID+II)))
              END DO
            END DO
          END DO
          IJBASE = IJBASE0 + ((ISTOP-KLA)/KLA+1) * ( KLA + JUMP )
          IIA = IIA0 + ((ISTOP-KLA)/KLA+1) * IINK
          IIB = IIB0 + ((ISTOP-KLA)/KLA+1) * IINK
          IIC = IIC0 - ((ISTOP-KLA)/KLA+1) * IINK
          IID = IID0 - ((ISTOP-KLA)/KLA+1) * IINK
!$acc end kernels
        END IF

        IF ( IIB <= IIC ) THEN
!$acc kernels
          IIBASE=0
          ZSIN45=SQRT(0.5)
!$acc loop independent
          DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
            DO IIL=0,KLA-1
              II = IIBASE + IJK * KINC3 + IIL
              IJ = IJBASE + IJK * KINC4 + IIL
              PC(IJA+IJ)=PA(IIA+II)+PA(IIB+II)
              PC(IJB+IJ)=ZSIN45*((PA(IIA+II)-PA(IIB+II))-(PB(IIA+II)+PB(IIB+II)))
              PC(IJC+IJ)=PB(IIB+II)-PB(IIA+II)
              PC(IJD+IJ)=-ZSIN45*((PA(IIA+II)-PA(IIB+II))+(PB(IIA+II)+PB(IIB+II)))
            END DO
          END DO
          IIBASE = IIBASE + KLA
          IJBASE = IJBASE + KLA
!$acc end kernels
        END IF

      ELSE ! KLA == M

!$acc kernels
!$acc loop independent
        DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
          DO IIL=0,KLA-1
            II = IIBASE + IJK * KINC3 + IIL
            IJ = IJBASE + IJK * KINC4 + IIL
            PC(IJA+IJ)=2.0*((PA(IIA+II)+PA(IIC+II))+PA(IIB+II))
            PC(IJB+IJ)=2.0*((PA(IIA+II)-PA(IIC+II))-PB(IIB+II))
            PC(IJC+IJ)=2.0*((PA(IIA+II)+PA(IIC+II))-PA(IIB+II))
            PC(IJD+IJ)=2.0*((PA(IIA+II)-PA(IIC+II))+PB(IIB+II))
          END DO
        END DO
        IIBASE = IIBASE + KLA
        IJBASE = IJBASE + KLA
!$acc end kernels
      END IF

      IBAD=0


      CASE ( 5 )
!
!     CODING FOR FACTOR 5
!     -------------------
      IIA=1
      IIB=IIA+2*IM-KLA
      IIC=IIB+2*IM
      IID=IIC
      IIE=IIB
      IJA=1
      IJB=IJA+JINK
      IJC=IJB+JINK
      IJD=IJC+JINK
      IJE=IJD+JINK

      IF ( KLA /=  IM ) THEN

!$acc kernels
!$acc loop independent
        DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
          DO IIL=0,KLA-1
            II = IIBASE + IJK * KINC3 + IIL
            IJ = IJBASE + IJK * KINC4 + IIL
            PC(IJA+IJ)=PA(IIA+II)+(PA(IIB+II)+PA(IIC+II))
            PC(IJB+IJ)=((PA(IIA+II)-0.25*(PA(IIB+II)+PA(IIC+II)))+XQRT5*(PA(IIB+II)-PA(IIC+ &
              II)))                                                        &
                -(XSIN72*PB(IIB+II)+XSIN36*PB(IIC+II))
            PC(IJC+IJ)=((PA(IIA+II)-0.25*(PA(IIB+II)+PA(IIC+II)))-XQRT5*(PA(IIB+II)-PA(IIC+ &
              II)))                                                        &
                -(XSIN36*PB(IIB+II)-XSIN72*PB(IIC+II))
            PC(IJD+IJ)=((PA(IIA+II)-0.25*(PA(IIB+II)+PA(IIC+II)))-XQRT5*(PA(IIB+II)-PA(IIC+ &
              II)))                                                        &
                +(XSIN36*PB(IIB+II)-XSIN72*PB(IIC+II))
            PC(IJE+IJ)=((PA(IIA+II)-0.25*(PA(IIB+II)+PA(IIC+II)))+XQRT5*(PA(IIB+II)-PA(IIC+ &
              II)))                                                        &
                +(XSIN72*PB(IIB+II)+XSIN36*PB(IIC+II))
          END DO
        END DO
        IIBASE = IIBASE + KLA
        IJBASE = IJBASE + KLA
!$acc end kernels
        IIA=IIA+IINK
        IINK=2*IINK
        IIB=IIB+IINK
        IIC=IIC+IINK
        IID=IID-IINK
        IIE=IIE-IINK
        IJBASE=IJBASE+JUMP
        JUMP=2*JUMP+JINK
        IF ( IIB /= IID ) THEN
!$acc kernels
          IJBASE0 = IJBASE
          IIA0 = IIA
          IIB0 = IIB
          IIC0 = IIC
          IID0 = IID
          IIE0 = IIE
!$acc loop independent private( IKB, IKC, IKD, IKE, ZC1, ZS1, ZC2, ZS2, ZC3, ZS3, ZC4, ZS4, IJBASE, IIA, IIB, IIC, IID, IIE )
          DO IIK=KLA,ISTOP,KLA
            IKB=IIK+IIK
            IKC=IKB+IKB
            IKD=IKC+IKB
            IKE=IKD+IKB
            ZC1=PTRIGS(IKB+1)
            ZS1=PTRIGS(IKB+2)
            ZC2=PTRIGS(IKC+1)
            ZS2=PTRIGS(IKC+2)
            ZC3=PTRIGS(IKD+1)
            ZS3=PTRIGS(IKD+2)
            ZC4=PTRIGS(IKE+1)
            ZS4=PTRIGS(IKE+2)
            IJBASE = IJBASE0 + (IIK-KLA)/KLA * (KLA + JUMP )
            IIA = IIA0 + (IIK-KLA)/KLA * IINK
            IIB = IIB0 + (IIK-KLA)/KLA * IINK
            IIC = IIC0 + (IIK-KLA)/KLA * IINK
            IID = IID0 - (IIK-KLA)/KLA * IINK
            IIE = IIE0 - (IIK-KLA)/KLA * IINK
!$acc loop independent
            DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ, ZA10, ZA11, ZA20, ZA21, ZB10, ZB11, ZB20, ZB21 )
              DO IIL=0,KLA-1
                II =          IJK * KINC3 + IIL
                IJ = IJBASE + IJK * KINC4 + IIL

                ZA10=(PA(IIA+II)-0.25*((PA(IIB+II)+PA(IIE+II))+(PA(IIC+II)+PA(IID+II)))) &
                    +XQRT5*((PA(IIB+II)+PA(IIE+II))-(PA(IIC+II)+PA(IID+II)))
                ZA20=(PA(IIA+II)-0.25*((PA(IIB+II)+PA(IIE+II))+(PA(IIC+II)+PA(IID+II)))) &
                    -XQRT5*((PA(IIB+II)+PA(IIE+II))-(PA(IIC+II)+PA(IID+II)))
                ZB10=(PB(IIA+II)-0.25*((PB(IIB+II)-PB(IIE+II))+(PB(IIC+II)-PB(IID+II)))) &
                    +XQRT5*((PB(IIB+II)-PB(IIE+II))-(PB(IIC+II)-PB(IID+II)))
                ZB20=(PB(IIA+II)-0.25*((PB(IIB+II)-PB(IIE+II))+(PB(IIC+II)-PB(IID+II)))) &
                    -XQRT5*((PB(IIB+II)-PB(IIE+II))-(PB(IIC+II)-PB(IID+II)))
                ZA11=XSIN72*(PB(IIB+II)+PB(IIE+II))+XSIN36*(PB(IIC+II)+PB(IID+II))
                ZA21=XSIN36*(PB(IIB+II)+PB(IIE+II))-XSIN72*(PB(IIC+II)+PB(IID+II))
                ZB11=XSIN72*(PA(IIB+II)-PA(IIE+II))+XSIN36*(PA(IIC+II)-PA(IID+II))
                ZB21=XSIN36*(PA(IIB+II)-PA(IIE+II))-XSIN72*(PA(IIC+II)-PA(IID+II))

                PC(IJA+IJ)=PA(IIA+II)+((PA(IIB+II)+PA(IIE+II))+(PA(IIC+II)+PA(IID+II)))
                PD(IJA+IJ)=PB(IIA+II)+((PB(IIB+II)-PB(IIE+II))+(PB(IIC+II)-PB(IID+II)))
                PC(IJB+IJ)=ZC1*(ZA10-ZA11)-ZS1*(ZB10+ZB11)
                PD(IJB+IJ)=ZS1*(ZA10-ZA11)+ZC1*(ZB10+ZB11)
                PC(IJE+IJ)=ZC4*(ZA10+ZA11)-ZS4*(ZB10-ZB11)
                PD(IJE+IJ)=ZS4*(ZA10+ZA11)+ZC4*(ZB10-ZB11)
                PC(IJC+IJ)=ZC2*(ZA20-ZA21)-ZS2*(ZB20+ZB21)
                PD(IJC+IJ)=ZS2*(ZA20-ZA21)+ZC2*(ZB20+ZB21)
                PC(IJD+IJ)=ZC3*(ZA20+ZA21)-ZS3*(ZB20-ZB21)
                PD(IJD+IJ)=ZS3*(ZA20+ZA21)+ZC3*(ZB20-ZB21)
              END DO
            END DO
          END DO
          IJBASE = IJBASE0 + ((ISTOP-KLA)/KLA+1) * ( KLA + JUMP )
          IIA = IIA0 + ((ISTOP-KLA)/KLA+1) * IINK
          IIB = IIB0 + ((ISTOP-KLA)/KLA+1) * IINK
          IIC = IIC0 + ((ISTOP-KLA)/KLA+1) * IINK
          IID = IID0 - ((ISTOP-KLA)/KLA+1) * IINK
          IIE = IIE0 - ((ISTOP-KLA)/KLA+1) * IINK
!$acc end kernels
        END IF

        IF ( IIB <= IID ) THEN
!$acc kernels
          IIBASE=0
!$acc loop independent
          DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
            DO IIL=0,KLA-1
              II = IIBASE + IJK * KINC3 + IIL
              IJ = IJBASE + IJK * KINC4 + IIL
              PC(IJA+IJ)=(PA(IIA+II)+PA(IIB+II))+PA(IIC+II)
              PC(IJB+IJ)=(XQRT5*(PA(IIA+II)-PA(IIB+II))+(0.25*(PA(IIA+II)+PA(IIB+II))-PA(IIC+ &
                II)))                                                        &
                  -(XSIN36*PB(IIA+II)+XSIN72*PB(IIB+II))
              PC(IJE+IJ)=-(XQRT5*(PA(IIA+II)-PA(IIB+II))+(0.25*(PA(IIA+II)+PA(IIB+II))-PA(IIC+ &
                II)))                                                         &
                  -(XSIN36*PB(IIA+II)+XSIN72*PB(IIB+II))
              PC(IJC+IJ)=(XQRT5*(PA(IIA+II)-PA(IIB+II))-(0.25*(PA(IIA+II)+PA(IIB+II))-PA(IIC+ &
                II)))                                                        &
                  -(XSIN72*PB(IIA+II)-XSIN36*PB(IIB+II))
              PC(IJD+IJ)=-(XQRT5*(PA(IIA+II)-PA(IIB+II))-(0.25*(PA(IIA+II)+PA(IIB+II))-PA(IIC+ &
                II)))                                                         &
                  -(XSIN72*PB(IIA+II)-XSIN36*PB(IIB+II))
            END DO
          END DO
          IIBASE = IIBASE + KLA
          IJBASE = IJBASE + KLA
!$acc end kernels
        END IF

      ELSE

!$acc kernels
        ZQQRT5=2.0*XQRT5
        ZSSIN36=2.0*XSIN36
        ZSSIN72=2.0*XSIN72
!$acc loop independent
        DO IJK=0,KLOT-1
!CDIR$ IVDEP
!!CDIR NODEP
!*VOCL LOOP,NOVREC
!$acc loop independent private( II, IJ )
          DO IIL=0,KLA-1
            II = IIBASE + IJK * KINC3 + IIL
            IJ = IJBASE + IJK * KINC4 + IIL
            PC(IJA+IJ)=2.0*(PA(IIA+II)+(PA(IIB+II)+PA(IIC+II)))
            PC(IJB+IJ)=(2.0*(PA(IIA+II)-0.25*(PA(IIB+II)+PA(IIC+II)))                 &
                +ZQQRT5*(PA(IIB+II)-PA(IIC+II)))-(ZSSIN72*PB(IIB+II)+ZSSIN36*PB(IIC+II))
            PC(IJC+IJ)=(2.0*(PA(IIA+II)-0.25*(PA(IIB+II)+PA(IIC+II)))                 &
                -ZQQRT5*(PA(IIB+II)-PA(IIC+II)))-(ZSSIN36*PB(IIB+II)-ZSSIN72*PB(IIC+II))
            PC(IJD+IJ)=(2.0*(PA(IIA+II)-0.25*(PA(IIB+II)+PA(IIC+II)))                 &
                -ZQQRT5*(PA(IIB+II)-PA(IIC+II)))+(ZSSIN36*PB(IIB+II)-ZSSIN72*PB(IIC+II))
            PC(IJE+IJ)=(2.0*(PA(IIA+II)-0.25*(PA(IIB+II)+PA(IIC+II)))                 &
                +ZQQRT5*(PA(IIB+II)-PA(IIC+II)))+(ZSSIN72*PB(IIB+II)+ZSSIN36*PB(IIC+II))
          END DO
        END DO
        IIBASE = IIBASE + KLA
        IJBASE = IJBASE + KLA