Skip to content
Snippets Groups Projects
fft.f90 70.5 KiB
Newer Older
  • Learn to ignore specific revisions
  • !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
    
      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)
    
      ! FIND FACTORS OF KN (8,6,5,4,3,2; ONLY ONE 8 ALLOWED)
    
      ! LOOK FOR SIXES FIRST, STORE FACTORS IN DESCENDING ORDER
    
        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
    
          IU=IU/IFAC
          IF ( IU == 1 ) EXIT
    
        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)
    
    SUBROUTINE FFT991( PA, PWORK, PTRIGS, PFAX, KJUMP, KN, KLOT, KSIGN )
    
      REAL,    DIMENSION(:),  INTENT(INOUT) :: PA
      REAL,    DIMENSION(:),  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
    
    !         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
    
    !     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 :: INB, IK, IFAC, IERR
      INTEGER :: IIBASE, IJBASE
    
      IF (MPPDB_INITIALIZED) THEN
        !Check all IN arrays
    
        CALL MPPDB_CHECK( PTRIGS, "FFT991 beg:PTRIGS" )
        CALL MPPDB_CHECK( PFAX, "FFT991 beg:KFAX" )
    
        CALL MPPDB_CHECK( PA, "FFT991 beg:PA" )
    
    !$acc data present( PA, PWORK )
    
      ! Initialisation of PWORK useful to compare results with MPPDB_CHECK (otherwise all values are not set by FFT991
      PWORK(:) = 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)
    
      IF(PFAX(10).NE.KN) CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'FFT991', 'PTRIGS and KFAX not initialised (set99 not yet called)' )
    
      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
    
          DO J=1,INVEX
            I = ISTART + ( J - 1 ) * KJUMP
            PA(I+1)=0.5*PA(I)
    
          IF ( MOD(KN,2) == 0 ) THEN
    
            DO J=1,INVEX
              I = I0 + ( J - 1 ) * KJUMP
              PA(I)=0.5*PA(I)
    
          DO IK=1,IFAX
            IFAC=PFAX(IK+1)
    
              CALL RPASSM( PA(IA:), PA(IA+ILA:), PWORK(:), PWORK(IFAC*ILA+1:), PTRIGS(:), &
                           KJUMP, INX, INVEX, KN, IFAC, ILA, IERR )
    
              CALL RPASSM( PWORK(:), PWORK(ILA+1:), PA(IA:), PA(IA+IFAC*ILA:), PTRIGS(:), &
                           INX, KJUMP, INVEX, KN, IFAC, ILA, IERR )
    
    !     IF NECESSARY, COPY RESULTS BACK TO PA
    !     -------------------------------------
          IF ( MOD(IFAX,2) == 1 ) THEN
    
    !$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) = PWORK(IIBASE:IIBASE+KN-1)
    
          PA(ISTART+KN::KJUMP) = 0.0
          PA(ISTART+KN+1::KJUMP) = 0.0
    
          ISTART=ISTART+INVEX*KJUMP
          INVEX=NVECLEN
    
      ELSE IF ( KSIGN == -1 ) THEN
    
    !     KSIGN=-1, GRIDPOINT TO SPECTRAL TRANSFORM
    
          DO IK=1,IFAX
            IFAC=PFAX(IFAX+2-IK)
    
              CALL QPASSM( PA(IA:), PA(IA+IFAC*ILA:), PWORK(:), PWORK(ILA+1:), PTRIGS(:), &
                           KJUMP, INX, INVEX, KN, IFAC, ILA, IERR )
    
              CALL QPASSM( PWORK(:), PWORK(IFAC*ILA+1:), PA(IA:), PA(IA+ILA:), PTRIGS(:), &
                           INX, KJUMP, INVEX, KN, IFAC, ILA, IERR )
    
    !     IF NECESSARY, COPY RESULTS BACK TO PA
    !     -------------------------------------
          IF ( MOD(IFAX,2) == 1 ) THEN
    
    !$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) = PWORK(IIBASE:IIBASE+KN-1)
    
    !     SHIFT PA(0) & FILL IN ZERO IMAG PARTS
    !     -------------------------------------
    
          DO J=1,INVEX
            IX  = ISTART + ( J - 1 ) * KJUMP
            PA(IX)=PA(IX+1)
            PA(IX+1)=0.0
    
          IF ( MOD(KN,2) == 0 ) THEN
    
            PA(ISTART+KN+1::KJUMP) = 0.0
    
          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')
    
    !$acc end data
    
      IF (MPPDB_INITIALIZED) THEN
        !Check all INOUT arrays
    
        CALL MPPDB_CHECK( PA, "FFT991 end:PA" )
    
        CALL MPPDB_CHECK( PWORK, "FFT991 end:PWORK" )
    
    SUBROUTINE RPASSM( PA, PB, PC, PD, PTRIGS, KINC3, KINC4, KLOT, KN, KFAC, KLA, KERR )
    
      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:
    
    !              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 )
    
          IM=KN/KFAC
          IINK=KLA
          JINK=KLA
          JUMP=(KFAC-1)*JINK
          ISTOP=(KN-KFAC)/(2*KFAC)
    
          IF (KLOT.GT.NVECLEN) GO TO 910
          IIBASE=0
          IJBASE=0
          IGO=KFAC-1
    
          SELECT CASE ( KFAC )
    
          IIA=1
          IIB=IIA+2*IM-KLA
          IJA=1
          IJB=IJA+JINK
    
          IF ( KLA /= IM ) THEN
    
    !$acc loop independent private( II, IJ )
              DO IJK=1,KLOT
                II = IIBASE + IIL - 1 + (IJK - 1 ) * KINC3
                IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                PC(IJA+IJ)=PA(IIA+II)+PA(IIB+II)
                PC(IJB+IJ)=PA(IIA+II)-PA(IIB+II)
    
            IIBASE=IIBASE+KLA
            IJBASE=IJBASE+KLA
    
            IIB=IIB-IINK
            IIBASE=0
            IJBASE=IJBASE+JUMP
    
            IF ( IIA /= IIB ) THEN
    
              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 private( II, IJ )
                  DO IJK=1,KLOT
                    II =         IIL - 1 + (IJK - 1 ) * KINC3
                    IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                    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))
    
              IJBASE = IJBASE0 + ((ISTOP-KLA)/KLA+1) * ( KLA + JUMP )
              IIA = IIA0 + ((ISTOP-KLA)/KLA+1) * IINK
              IIB = IIB0 - ((ISTOP-KLA)/KLA+1) * IINK
    
            IF ( IIA <= IIB ) THEN
    
    !$acc loop independent private( II, IJ )
                DO IJK=1,KLOT
                  II = IIBASE + IIL - 1 + (IJK - 1 ) * KINC3
                  IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                  PC(IJA+IJ)=PA(IIA+II)
                  PC(IJB+IJ)=-PB(IIA+II)
    
              IIBASE=IIBASE+KLA
              IJBASE=IJBASE+KLA
    
    !$acc loop independent private( II, IJ )
              DO IJK=1,KLOT
                II = IIBASE + IIL - 1 + (IJK - 1 ) * KINC3
                IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                PC(IJA+IJ)=2.0*(PA(IIA+II)+PA(IIB+II))
                PC(IJB+IJ)=2.0*(PA(IIA+II)-PA(IIB+II))
    
            IIBASE=IIBASE+KLA
            IJBASE=IJBASE+KLA
    
          IIA=1
          IIB=IIA+2*IM-KLA
          IIC=IIB
          IJA=1
          IJB=IJA+JINK
          IJC=IJB+JINK
    
          IF ( KLA /= IM ) THEN
    
    !$acc loop independent private( II, IJ )
              DO IJK=1,KLOT
                II = IIBASE + IIL - 1 + (IJK - 1 ) * KINC3
                IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                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)))
    
            IIBASE=IIBASE+KLA
            IJBASE=IJBASE+KLA
    
            IIB=IIB+IINK
            IIC=IIC-IINK
            IJBASE=IJBASE+JUMP
    
            IF ( IIA /= IIC ) THEN
    
              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 private( II, IJ )
                  DO IJK=1,KLOT
                    II =         IIL - 1 + (IJK - 1 ) * KINC3
                    IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                    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))))
    
              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
    
            IF ( IIA <= IIC ) THEN
    
    !$acc loop independent private( II, IJ )
                DO IJK=1,KLOT
                  II = IIBASE + IIL - 1 + (IJK - 1 ) * KINC3
                  IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                  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))
    
              IIBASE = IIBASE + KLA
              IJBASE = IJBASE + KLA
    
    !$acc loop independent private( II, IJ )
              DO IJK=1,KLOT
                II = IIBASE + IIL - 1 + (IJK - 1 ) * KINC3
                IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                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))
    
            IIBASE = IIBASE + KLA
            IJBASE = IJBASE + KLA
    
          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 loop independent private( II, IJ )
              DO IJK=1,KLOT
                II = IIBASE + IIL - 1 + (IJK - 1 ) * KINC3
                IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                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)
    
            IIBASE = IIBASE + KLA
            IJBASE = IJBASE + KLA
    
            IIB=IIB+IINK
            IIC=IIC-IINK
            IID=IID-IINK
            IJBASE=IJBASE+JUMP
    
            IF ( IIB /=  IIC ) THEN
    
              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 private( II, IJ )
                  DO IJK=1,KLOT
                    II =         IIL - 1 + (IJK - 1 ) * KINC3
                    IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                    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)))
    
              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
    
            IF ( IIB <= IIC ) THEN
    
              IIBASE=0
              ZSIN45=SQRT(0.5)
    
    !$acc loop independent private( II, IJ )
                DO IJK=1,KLOT
                  II = IIBASE + IIL - 1 + (IJK - 1 ) * KINC3
                  IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                  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)))
    
              IIBASE = IIBASE + KLA
              IJBASE = IJBASE + KLA
    
    !$acc loop independent private( II, IJ )
              DO IJK=1,KLOT
                II = IIBASE + IIL - 1 + (IJK - 1 ) * KINC3
                IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                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))
    
            IIBASE = IIBASE + KLA
            IJBASE = IJBASE + KLA
    
          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 loop independent private( II, IJ )
              DO IJK=1,KLOT
                II = IIBASE + IIL - 1 + (IJK - 1 ) * KINC3
                IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                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))
    
            IIBASE = IIBASE + KLA
            IJBASE = IJBASE + KLA
    
            IIB=IIB+IINK
            IIC=IIC+IINK
            IID=IID-IINK
            IIE=IIE-IINK
            IJBASE=IJBASE+JUMP
    
            IF ( IIB /= IID ) THEN
    
              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 private( II, IJ, ZA10, ZA11, ZA20, ZA21, ZB10, ZB11, ZB20, ZB21 )
                  DO IJK=1,KLOT
                    II =         IIL - 1 + (IJK - 1 ) * KINC3
                    IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
    
                    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)
    
              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
    
            IF ( IIB <= IID ) THEN
    
    !$acc loop independent private( II, IJ )
                DO IJK=1,KLOT
                  II = IIBASE + IIL - 1 + (IJK - 1 ) * KINC3
                  IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                  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))
    
              IIBASE = IIBASE + KLA
              IJBASE = IJBASE + KLA
    
            ZQQRT5=2.0*XQRT5
            ZSSIN36=2.0*XSIN36
            ZSSIN72=2.0*XSIN72
    
    !$acc loop independent private( II, IJ )
              DO IJK=1,KLOT
                II = IIBASE + IIL - 1 + (IJK - 1 ) * KINC3
                IJ = IJBASE + IIL - 1 + (IJK - 1 ) * KINC4
                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))
    
            IIBASE = IIBASE + KLA
            IJBASE = IJBASE + KLA