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, KSZA, KSZW, KSZT )
  USE MODE_MPPDB

  IMPLICIT NONE

      REAL,    DIMENSION(KSZA), INTENT(INOUT) :: PA
      REAL,    DIMENSION(KSZW), INTENT(OUT)   :: PWORK
      REAL,    DIMENSION(KSZT), INTENT(IN)    :: PTRIGS
      INTEGER, DIMENSION(19),   INTENT(IN)    :: PFAX
      INTEGER,                  INTENT(IN)    :: KJUMP, KN, KLOT, KSIGN
      INTEGER,                  INTENT(IN)    :: KSZA, KSZW, KSZT
!
!     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

  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

!$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.

#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