MODULE MDL
IMPLICIT NONE
INTEGER,PARAMETER :: RK=SELECTED_REAL_KIND(15,307)
INTEGER,PARAMETER :: NUM=2**10
INTEGER,PARAMETER :: FLA=1
REAL(RK),PARAMETER :: PI=3.141592653589793238460_RK
REAL(RK),PARAMETER :: TWO_PI=2.0_RK*PI
CONTAINS
! DISCRETE FOURIER TRANSFORM (NRF77)
PURE SUBROUTINE FFT_(NUM,DIR,ARR)
!$ACC ROUTINE SEQ
INTEGER,INTENT(IN) :: NUM
INTEGER,INTENT(IN) :: DIR
REAL(RK),DIMENSION(2*NUM),INTENT(INOUT) :: ARR
INTEGER :: N,I,J,M,LIM,STE
REAL(RK) :: PIM,PRE,ANG,WR,WI,WPR,WPI,WX,MUL
N=2*NUM
J=1
MUL = REAL(DIR,RK)*TWO_PI
DO I=1,N,2
IF(J.GT.I)THEN
PRE=ARR(J)
PIM=ARR(J+1)
ARR(J)=ARR(I)
ARR(J+1)=ARR(I+1)
ARR(I)=PRE
ARR(I+1)=PIM
ENDIF
M=N/2
1 IF ((M.GE.2).AND.(J.GT.M)) THEN
J=J-M
M=M/2
GOTO 1
ENDIF
J=J+M
ENDDO
LIM=2
2 IF (N.GT.LIM) THEN
STE=2*LIM
ANG=MUL/REAL(LIM,RK)
WPI=SIN(ANG)
ANG=SIN(0.5_RK*ANG)
WPR=-2.0_RK*ANG*ANG
WR=1.0_RK
WI=0.0_RK
DO M=1,LIM,2
DO I=M,N,STE
J=I+LIM
PRE=WR*ARR(J)-WI*ARR(J+1)
PIM=WR*ARR(J+1)+WI*ARR(J)
ARR(J)=ARR(I)-PRE
ARR(J+1)=ARR(I+1)-PIM
ARR(I)=ARR(I)+PRE
ARR(I+1)=ARR(I+1)+PIM
ENDDO
WX=WR
WR=WR*WPR-WI*WPI+WR
WI=WI*WPR+WX*WPI+WI
ENDDO
LIM=STE
GOTO 2
ENDIF
END SUBROUTINE FFT_
! DISCRETE (LINEAR) FRACTIONAL FOURIER TRANSFORM
PURE SUBROUTINE FFRFT_(NUM,PAR,ARR)
!$ACC ROUTINE SEQ
INTEGER,INTENT(IN) :: NUM
REAL(RK),INTENT(IN) :: PAR
REAL(RK),DIMENSION(2*NUM),INTENT(INOUT) :: ARR
INTEGER :: I
REAL(RK) :: FAC
REAL(RK),DIMENSION(NUM) :: MUL,COS_MUL,SIN_MUL
REAL(RK),DIMENSION(4*NUM) :: ONE,TWO,TRE
REAL(RK),DIMENSION(2*NUM) :: TMP
FAC=PAR*PI/REAL(NUM,RK)
MUL=FAC*REAL([(I-1,I=1,NUM,1)],RK)**2
COS_MUL=COS(MUL)
SIN_MUL=SIN(MUL)
ONE(1:2*NUM:2)=ARR(1:2*NUM:2)*COS_MUL-ARR(2:2*NUM:2)*SIN_MUL
ONE(2:2*NUM:2)=ARR(1:2*NUM:2)*SIN_MUL+ARR(2:2*NUM:2)*COS_MUL
TWO(1:2*NUM:2) = +COS_MUL
TWO(2:2*NUM:2) = -SIN_MUL
ONE(2*NUM+1:4*NUM:1)=0.0_RK
MUL=-FAC*REAL([(I-1-2*NUM,I=NUM+1,2*NUM,1)],RK)**2
TWO(2*NUM+1:4*NUM:2)=COS(MUL)
TWO(2*NUM+2:4*NUM:2)=SIN(MUL)
CALL FFT_(2*NUM,+1,ONE)
CALL FFT_(2*NUM,+1,TWO)
TRE=ONE
ONE(1:4*NUM:2)=TRE(1:4*NUM:2)*TWO(1:4*NUM:2)-TRE(2:4*NUM:2)*TWO(2:4*NUM:2)
ONE(2:4*NUM:2)=TRE(1:4*NUM:2)*TWO(2:4*NUM:2)+TRE(2:4*NUM:2)*TWO(1:4*NUM:2)
CALL FFT_(2*NUM,-1,ONE)
ARR=ONE(1:2*NUM:1)/REAL(2*NUM,RK)
TMP=ARR
ARR(1:2*NUM:2)=TMP(1:2*NUM:2)*COS_MUL-TMP(2:2*NUM:2)*SIN_MUL
ARR(2:2*NUM:2)=TMP(1:2*NUM:2)*SIN_MUL+TMP(2:2*NUM:2)*COS_MUL
END SUBROUTINE FFRFT_
END MODULE MDL