tpm_fftw.F90 Source File


Files dependent on this one

sourcefile~~tpm_fftw.f90~~AfferentGraph sourcefile~tpm_fftw.f90 tpm_fftw.F90 sourcefile~dealloc_resol_mod.f90~2 dealloc_resol_mod.F90 sourcefile~dealloc_resol_mod.f90~2->sourcefile~tpm_fftw.f90 sourcefile~ftdir_mod.f90~2 ftdir_mod.F90 sourcefile~ftdir_mod.f90~2->sourcefile~tpm_fftw.f90 sourcefile~ftdirad_mod.f90 ftdirad_mod.F90 sourcefile~ftdirad_mod.f90->sourcefile~tpm_fftw.f90 sourcefile~ftinv_mod.f90~2 ftinv_mod.F90 sourcefile~ftinv_mod.f90~2->sourcefile~tpm_fftw.f90 sourcefile~ftinvad_mod.f90 ftinvad_mod.F90 sourcefile~ftinvad_mod.f90->sourcefile~tpm_fftw.f90 sourcefile~set_resol_mod.f90~2 set_resol_mod.F90 sourcefile~set_resol_mod.f90~2->sourcefile~tpm_fftw.f90 sourcefile~setup_trans.f90 setup_trans.F90 sourcefile~setup_trans.f90->sourcefile~tpm_fftw.f90 sourcefile~setup_trans.f90~2 setup_trans.F90 sourcefile~setup_trans.f90~2->sourcefile~tpm_fftw.f90 sourcefile~trans_end.f90~2 trans_end.F90 sourcefile~trans_end.f90~2->sourcefile~tpm_fftw.f90 sourcefile~ftdir_ctlad_mod.f90 ftdir_ctlad_mod.F90 sourcefile~ftdir_ctlad_mod.f90->sourcefile~ftdirad_mod.f90 sourcefile~ftinv_ctlad_mod.f90 ftinv_ctlad_mod.F90 sourcefile~ftinv_ctlad_mod.f90->sourcefile~ftinvad_mod.f90 sourcefile~dir_trans_ctlad_mod.f90 dir_trans_ctlad_mod.F90 sourcefile~dir_trans_ctlad_mod.f90->sourcefile~ftdir_ctlad_mod.f90 sourcefile~inv_trans_ctlad_mod.f90 inv_trans_ctlad_mod.F90 sourcefile~inv_trans_ctlad_mod.f90->sourcefile~ftinv_ctlad_mod.f90 sourcefile~dir_transad.f90~2 dir_transad.F90 sourcefile~dir_transad.f90~2->sourcefile~dir_trans_ctlad_mod.f90 sourcefile~inv_transad.f90~2 inv_transad.F90 sourcefile~inv_transad.f90~2->sourcefile~inv_trans_ctlad_mod.f90

Source Code

! (C) Copyright 2000- ECMWF.
! (C) Copyright 2000- Meteo-France.
! 
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!

MODULE TPM_FFTW
!   Author.
!   -------
!     George Mozdzynski
!
!   Modifications.
!   -------------- 
!     Original      October 2014
!     R. El Khatib 01-Sep-2015 More subroutines for better modularity

USE, INTRINSIC :: ISO_C_BINDING

USE PARKIND1   ,ONLY : JPIB, JPIM, JPRB, JPRD
USE MPL_MODULE ,ONLY : MPL_MYRANK
USE YOMHOOK    ,ONLY : LHOOK, DR_HOOK, JPHOOK

IMPLICIT NONE

SAVE

#ifdef __NEC__
! From NLC (NEC Numeric Library Collection)
#include "aslfftw3.f03"
#define FFTW_NO_SIMD 0
#else
#include "fftw3.f03"
#endif

PRIVATE
PUBLIC CREATE_PLAN_FFTW, DESTROY_PLAN_FFTW, DESTROY_PLANS_FFTW, INIT_PLANS_FFTW, &
      & FFTW_RESOL, TW, EXEC_FFTW, EXEC_EFFTW

TYPE FFTW_TYPE
  INTEGER(KIND=JPIM),ALLOCATABLE :: N_PLANS(:)
  TYPE(FFTW_PLAN),POINTER :: FFTW_PLANS(:) => NULL()
  INTEGER(KIND=JPIM) :: N_MAX=0         ! maximum number of latitudes
  INTEGER(KIND=JPIM) :: N_MAX_PLANS=4   ! maximum number of plans for each active latitudes
END TYPE FFTW_TYPE


TYPE FFTW_PLAN
  INTEGER(KIND=JPIM) :: NPLAN_ID=123456
  INTEGER(KIND=JPIB) :: NPLAN
  INTEGER(KIND=JPIM) :: NLOT
  INTEGER(KIND=JPIM) :: NTYPE
  TYPE(FFTW_PLAN),POINTER :: NEXT_PLAN => NULL()
END TYPE FFTW_PLAN

TYPE(FFTW_TYPE),ALLOCATABLE,TARGET :: FFTW_RESOL(:)
TYPE(FFTW_TYPE),POINTER     :: TW



! ------------------------------------------------------------------
CONTAINS
! ------------------------------------------------------------------


SUBROUTINE INIT_PLANS_FFTW(KDLON)
INTEGER(KIND=JPIM),INTENT(IN) :: KDLON

#include "abor1.intfb.h"

TW%N_MAX=KDLON
ALLOCATE(TW%FFTW_PLANS(TW%N_MAX))
ALLOCATE(TW%N_PLANS(TW%N_MAX))
TW%N_PLANS(:)=0
RETURN  
END SUBROUTINE INIT_PLANS_FFTW


SUBROUTINE CREATE_PLAN_FFTW(KPLAN,KTYPE,KN,KLOT)
INTEGER(KIND=JPIB),INTENT(OUT) :: KPLAN
INTEGER(KIND=JPIM),INTENT(IN) :: KTYPE,KN,KLOT

INTEGER(KIND=JPIB) :: IPLAN
INTEGER(KIND=JPIM) :: IRANK, ISTRIDE
INTEGER(KIND=JPIM) :: JL
INTEGER(KIND=JPIM) :: IRDIST,ICDIST,IN(1),IEMBED(1)
REAL(KIND=JPRB), POINTER :: ZDUM(:)
TYPE(C_PTR) :: ZDUMP
LOGICAL :: LLFOUND
LOGICAL :: LLRESTRICT_PLANS=.TRUE.
TYPE(FFTW_PLAN),POINTER :: CURR_FFTW_PLAN,START_FFTW_PLAN
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE, ZHOOK_HANDLE2
IF (LHOOK) CALL DR_HOOK('CREATE_PLAN_FFTW',0,ZHOOK_HANDLE)

IF( KN > TW%N_MAX )THEN
  CALL ABOR1('CREATE_PLAN_FFTW: KN > N_MAX THAT WAS INITIALISED IN INIT_PLANS_FFTW')
ENDIF

IRANK=1
ISTRIDE=1
IN(1)=KN
IEMBED(1)=IN(1)
ICDIST=KN/2+1
IRDIST=ICDIST*2

!$OMP CRITICAL (FFTW_CREATE)
LLFOUND=.FALSE.
IF( TW%FFTW_PLANS(KN)%NPLAN_ID /= 123456 )THEN
  WRITE(*,'("CREATE_PLAN_FFTW.1: PLAN_ID=",I10)')TW%FFTW_PLANS(KN)%NPLAN_ID
  CALL ABOR1('CREATE_PLAN_FFTW.1: NPLAN_ID /= 123456')
ENDIF
CURR_FFTW_PLAN=>TW%FFTW_PLANS(KN)
IF( CURR_FFTW_PLAN%NPLAN_ID /= 123456 )THEN
  WRITE(*,'("CREATE_PLAN_FFTW.2: PLAN_ID=",I10)')CURR_FFTW_PLAN%NPLAN_ID
  CALL ABOR1('CREATE_PLAN_FFTW.2: NPLAN_ID /= 123456')
ENDIF
! search for plan in existing plans
DO JL=1,TW%N_PLANS(KN)
  IF( KLOT == CURR_FFTW_PLAN%NLOT .AND. KTYPE == CURR_FFTW_PLAN%NTYPE )THEN
    LLFOUND=.TRUE.
    IPLAN=CURR_FFTW_PLAN%NPLAN
    EXIT
  ELSEIF( JL /= TW%N_PLANS(KN) )THEN
    CURR_FFTW_PLAN=>CURR_FFTW_PLAN%NEXT_PLAN
    IF( CURR_FFTW_PLAN%NPLAN_ID /= 123456 )THEN
      WRITE(*,'("CREATE_PLAN_FFTW.3: PLAN_ID=",I10)')CURR_FFTW_PLAN%NPLAN_ID
      CALL ABOR1('CREATE_PLAN_FFTW.3: NPLAN_ID /= 123456')
    ENDIF
  ENDIF
ENDDO
IF( .NOT.LLFOUND )THEN
  IF( LLRESTRICT_PLANS )THEN
    IF( TW%N_PLANS(KN) == TW%N_MAX_PLANS )THEN
      ! destroy the plan at the start of the list
!     WRITE(*,'("CREATE_PLAN_FFTW: BEG: DESTROYING A PLAN AT THE START OF THE LIST")')
       IF (JPRB == JPRD) THEN
          CALL DFFTW_DESTROY_PLAN(TW%FFTW_PLANS(KN)%NPLAN)
       ELSE
          CALL SFFTW_DESTROY_PLAN(TW%FFTW_PLANS(KN)%NPLAN)
       END IF
      TW%FFTW_PLANS(KN)%NPLAN_ID=999999
      START_FFTW_PLAN=>TW%FFTW_PLANS(KN)
      TW%FFTW_PLANS(KN)=TW%FFTW_PLANS(KN)%NEXT_PLAN
      ! DEALLOCATE(START_FFTW_PLAN)
      TW%N_PLANS(KN)=TW%N_PLANS(KN)-1
!     WRITE(*,'("CREATE_PLAN_FFTW: END: DESTROYING A PLAN AT THE START OF THE LIST")')
    ENDIF
  ENDIF
  IF (JPRB == JPRD) THEN
     ZDUMP=FFTW_ALLOC_COMPLEX(INT(1,C_SIZE_T))
  ELSE
     ZDUMP=FFTWF_ALLOC_COMPLEX(INT(1,C_SIZE_T))
  END IF
  CALL C_F_POINTER(ZDUMP,ZDUM,[2])
  IF( KTYPE==1 )THEN
     IF (LHOOK) CALL DR_HOOK('FFTW_PLAN_MANY_DFT_C2R',0,ZHOOK_HANDLE2)
     IF (JPRB == JPRD) THEN
        CALL DFFTW_PLAN_MANY_DFT_C2R(IPLAN,IRANK,IN,KLOT,ZDUM,IEMBED,ISTRIDE,ICDIST,&
             & ZDUM,IEMBED,ISTRIDE,IRDIST,FFTW_ESTIMATE+FFTW_NO_SIMD)
     ELSE
        CALL SFFTW_PLAN_MANY_DFT_C2R(IPLAN,IRANK,IN,KLOT,ZDUM,IEMBED,ISTRIDE,ICDIST,&
             & ZDUM,IEMBED,ISTRIDE,IRDIST,FFTW_ESTIMATE+FFTW_NO_SIMD)
     END IF
     IF (LHOOK) CALL DR_HOOK('FFTW_PLAN_MANY_DFT_C2R',1,ZHOOK_HANDLE2)
  ELSEIF( KTYPE==-1 )THEN
     IF (LHOOK) CALL DR_HOOK('FFTW_PLAN_MANY_DFT_R2C',0,ZHOOK_HANDLE2)
     IF (JPRB == JPRD) THEN
        CALL DFFTW_PLAN_MANY_DFT_R2C(IPLAN,IRANK,IN,KLOT,ZDUM,IEMBED,ISTRIDE,IRDIST,&
             & ZDUM,IEMBED,ISTRIDE,ICDIST,FFTW_ESTIMATE+FFTW_NO_SIMD)
     ELSE
        CALL SFFTW_PLAN_MANY_DFT_R2C(IPLAN,IRANK,IN,KLOT,ZDUM,IEMBED,ISTRIDE,IRDIST,&
             & ZDUM,IEMBED,ISTRIDE,ICDIST,FFTW_ESTIMATE+FFTW_NO_SIMD)       
     END IF
     IF (LHOOK) CALL DR_HOOK('FFTW_PLAN_MANY_DFT_R2C',1,ZHOOK_HANDLE2)
  ELSE
    CALL ABOR1('FFTW_PLAN: INVALID KTYPE')
  ENDIF
  IF (JPRB == JPRD) THEN
     CALL FFTW_FREE(ZDUMP)
  ELSE
     CALL FFTWF_FREE(ZDUMP)
  END IF
  KPLAN=IPLAN
  TW%N_PLANS(KN)=TW%N_PLANS(KN)+1
  IF( TW%N_PLANS(KN) /= 1 )THEN
    ALLOCATE(CURR_FFTW_PLAN%NEXT_PLAN)
    CURR_FFTW_PLAN=>CURR_FFTW_PLAN%NEXT_PLAN
  ENDIF
  IF( CURR_FFTW_PLAN%NPLAN_ID /= 123456 )THEN
    WRITE(*,'("CREATE_PLAN_FFTW.4: PLAN_ID=",I10)')CURR_FFTW_PLAN%NPLAN_ID
    CALL ABOR1('CREATE_PLAN_FFTW.4: NPLAN_ID /= 123456')
  ENDIF
  CURR_FFTW_PLAN%NPLAN=IPLAN
  CURR_FFTW_PLAN%NLOT=KLOT
  CURR_FFTW_PLAN%NTYPE=KTYPE
  CURR_FFTW_PLAN%NEXT_PLAN=>NULL()
! write(*,'("CREATE_PLAN_FFTW: KN=",I5," NPLANS=",I3," KLOT=",I6," KTYPE=",I2,&
!  & " NEW IPLAN=",Z16)')KN,TW%N_PLANS(KN),KLOT,KTYPE,IPLAN
ELSE
  KPLAN=IPLAN
! write(*,'("CREATE_PLAN_FFTW: KN=",I5," NPLANS=",I3," KLOT=",I6," KTYPE=",I2,&
!  & " CUR IPLAN=",Z16)')KN,TW%N_PLANS(KN),KLOT,KTYPE,IPLAN
ENDIF
!$OMP END CRITICAL (FFTW_CREATE)

IF (LHOOK) CALL DR_HOOK('CREATE_PLAN_FFTW',1,ZHOOK_HANDLE)
RETURN
END SUBROUTINE CREATE_PLAN_FFTW


SUBROUTINE DESTROY_PLAN_FFTW(KPLAN)
INTEGER(KIND=JPIB),INTENT(IN) :: KPLAN
!$OMP CRITICAL (FFTW_DESTROY)
IF (JPRB == JPRD) THEN
   CALL DFFTW_DESTROY_PLAN(KPLAN)
ELSE
   CALL SFFTW_DESTROY_PLAN(KPLAN)
END IF
!$OMP END CRITICAL (FFTW_DESTROY)
RETURN
END SUBROUTINE DESTROY_PLAN_FFTW


SUBROUTINE DESTROY_PLANS_FFTW
INTEGER(KIND=JPIM) :: JL, JN
TYPE(FFTW_PLAN),POINTER :: CURR_FFTW_PLAN, NEXT_FFTW_PLAN
DO JN=1,TW%N_MAX
  CURR_FFTW_PLAN=>TW%FFTW_PLANS(JN)
  DO JL=1,TW%N_PLANS(JN)
    CALL DESTROY_PLAN_FFTW(CURR_FFTW_PLAN%NPLAN)
    NEXT_FFTW_PLAN=>CURR_FFTW_PLAN%NEXT_PLAN
    IF( JL /= 1 ) THEN
      DEALLOCATE( CURR_FFTW_PLAN )
    ENDIF
    CURR_FFTW_PLAN => NEXT_FFTW_PLAN
  ENDDO
ENDDO
IF( ASSOCIATED(TW) ) THEN
  IF( ASSOCIATED(TW%FFTW_PLANS) )  DEALLOCATE(TW%FFTW_PLANS)
  IF( ALLOCATED(TW%N_PLANS) )     DEALLOCATE(TW%N_PLANS)
  TW%N_MAX=0
ENDIF
RETURN
END SUBROUTINE DESTROY_PLANS_FFTW

SUBROUTINE EXEC_FFTW(KTYPE,KRLEN,KCLEN,KOFF,KFIELDS,LD_ALL,PREEL)

INTEGER(KIND=JPIM),INTENT(IN)   :: KTYPE
INTEGER(KIND=JPIM),INTENT(IN)   :: KRLEN
INTEGER(KIND=JPIM),INTENT(IN)   :: KCLEN
INTEGER(KIND=JPIM),INTENT(IN)   :: KOFF
INTEGER(KIND=JPIM),INTENT(IN)   :: KFIELDS
LOGICAL           ,INTENT(IN)   :: LD_ALL
REAL(KIND=JPRB), INTENT(INOUT)  :: PREEL(:,:)

REAL(KIND=JPRB), POINTER :: ZFFT(:,:)
REAL(KIND=JPRB), POINTER :: ZFFT1(:)
TYPE(C_PTR) :: ZFFTP, ZFFT1P

INTEGER(KIND=JPIM) :: JJ,JF

INTEGER(KIND=JPIB) :: IPLAN_C2R, IPLAN_C2R1
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE, ZHOOK_HANDLE2

#include "abor1.intfb.h"

IF (LHOOK) CALL DR_HOOK('EXEC_FFTW',0,ZHOOK_HANDLE)

IF ( (KTYPE /= -1) .AND. (KTYPE /=1) ) THEN
  CALL ABOR1('TPM_FFTW:EXEC_FFTW : WRONG VALUE KTYPE')
ENDIF

IF( LD_ALL )THEN
  CALL CREATE_PLAN_FFTW(IPLAN_C2R,KTYPE,KRLEN,KFIELDS)
  IF (JPRB == JPRD) THEN
     ZFFTP=FFTW_ALLOC_COMPLEX(INT(KCLEN/2*KFIELDS,C_SIZE_T))
  ELSE
     ZFFTP=FFTWF_ALLOC_COMPLEX(INT(KCLEN/2*KFIELDS,C_SIZE_T))
  END IF
  CALL C_F_POINTER(ZFFTP,ZFFT,[KCLEN,KFIELDS])
  IF (KTYPE==1) THEN
    DO JF=1,KFIELDS
      DO JJ=1,KCLEN
        ZFFT(JJ,JF) =PREEL(JF,KOFF+JJ-1)
      ENDDO
    ENDDO
    IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_C2R',0,ZHOOK_HANDLE2)
    IF (JPRB == JPRD) THEN
       CALL DFFTW_EXECUTE_DFT_C2R(IPLAN_C2R,ZFFT,ZFFT)
    ELSE
       CALL SFFTW_EXECUTE_DFT_C2R(IPLAN_C2R,ZFFT,ZFFT)
    END IF
    IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_C2R',1,ZHOOK_HANDLE2)
    DO JJ=1,KRLEN
      DO JF=1,KFIELDS
        PREEL(JF,KOFF+JJ-1)=ZFFT(JJ,JF)
      ENDDO
    ENDDO
  ELSE
    DO JF=1,KFIELDS
      DO JJ=1,KRLEN
        ZFFT(JJ,JF) =PREEL(JF,KOFF+JJ-1)
      ENDDO
    ENDDO
    IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_R2C',0,ZHOOK_HANDLE2)
    IF (JPRB == JPRD) THEN
       CALL DFFTW_EXECUTE_DFT_R2C(IPLAN_C2R,ZFFT,ZFFT)
    ELSE
       CALL SFFTW_EXECUTE_DFT_R2C(IPLAN_C2R,ZFFT,ZFFT)
    END IF
    IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_R2C',1,ZHOOK_HANDLE2)
    DO JJ=1,KCLEN
      DO JF=1,KFIELDS
        PREEL(JF,KOFF+JJ-1)=ZFFT(JJ,JF)/REAL(KRLEN,JPRB)
      ENDDO
    ENDDO
  ENDIF
  IF (JPRB == JPRD) THEN
     CALL FFTW_FREE(ZFFTP)
  ELSE
     CALL FFTWF_FREE(ZFFTP)
  END IF
ELSE
  CALL CREATE_PLAN_FFTW(IPLAN_C2R1,KTYPE,KRLEN,1)
  IF (JPRB == JPRD) THEN
     ZFFT1P=FFTW_ALLOC_COMPLEX(INT(KCLEN/2,C_SIZE_T))
  ELSE
     ZFFT1P=FFTWF_ALLOC_COMPLEX(INT(KCLEN/2,C_SIZE_T))
  END IF
  CALL C_F_POINTER(ZFFT1P,ZFFT1,[KCLEN])
  IF (KTYPE==1) THEN
    DO JF=1,KFIELDS
      DO JJ=1,KCLEN
        ZFFT1(JJ) =PREEL(JF,KOFF+JJ-1)
      ENDDO
      IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_C2R',0,ZHOOK_HANDLE2)
      IF (JPRB == JPRD) THEN
         CALL DFFTW_EXECUTE_DFT_C2R(IPLAN_C2R1,ZFFT1,ZFFT1)
      ELSE
         CALL SFFTW_EXECUTE_DFT_C2R(IPLAN_C2R1,ZFFT1,ZFFT1)
      END IF
      IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_C2R',1,ZHOOK_HANDLE2)
      DO JJ=1,KRLEN
        PREEL(JF,KOFF+JJ-1)=ZFFT1(JJ)
      ENDDO
    ENDDO
  ELSE
    DO JF=1,KFIELDS
      DO JJ=1,KRLEN
        ZFFT1(JJ) =PREEL(JF,KOFF+JJ-1)
      ENDDO
      IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_R2C',0,ZHOOK_HANDLE2)
      IF (JPRB == JPRD) THEN
         CALL DFFTW_EXECUTE_DFT_R2C(IPLAN_C2R1,ZFFT1,ZFFT1)
      ELSE
         CALL SFFTW_EXECUTE_DFT_R2C(IPLAN_C2R1,ZFFT1,ZFFT1)
      END IF
      IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_R2C',1,ZHOOK_HANDLE2)
      DO JJ=1,KCLEN
        PREEL(JF,KOFF+JJ-1)=ZFFT1(JJ)/REAL(KRLEN,JPRB)
      ENDDO
    ENDDO
  ENDIF
  IF (JPRB == JPRD) THEN
     CALL FFTW_FREE(ZFFT1P)
  ELSE 
     CALL FFTWF_FREE(ZFFT1P)
  END IF  
ENDIF

IF (LHOOK) CALL DR_HOOK('EXEC_FFTW',1,ZHOOK_HANDLE)
END SUBROUTINE EXEC_FFTW

SUBROUTINE EXEC_EFFTW(KTYPE,KRLEN,KCLEN,KOFF,KFIELDS,LD_ALL,PREEL)

INTEGER(KIND=JPIM),INTENT(IN)   :: KTYPE
INTEGER(KIND=JPIM),INTENT(IN)   :: KRLEN
INTEGER(KIND=JPIM),INTENT(IN)   :: KCLEN
INTEGER(KIND=JPIM),INTENT(IN)   :: KOFF
INTEGER(KIND=JPIM),INTENT(IN)   :: KFIELDS
LOGICAL           ,INTENT(IN)   :: LD_ALL
REAL(KIND=JPRB), INTENT(INOUT)  :: PREEL(:,:)

REAL(KIND=JPRB), POINTER :: ZFFT(:,:)
REAL(KIND=JPRB), POINTER :: ZFFT1(:)
TYPE(C_PTR) :: ZFFTP, ZFFT1P

INTEGER(KIND=JPIM) :: JJ,JF

INTEGER(KIND=JPIB) :: IPLAN_C2R, IPLAN_C2R1
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE, ZHOOK_HANDLE2

#include "abor1.intfb.h"

IF (LHOOK) CALL DR_HOOK('EXEC_EFFTW',0,ZHOOK_HANDLE)

IF ( (KTYPE /= -1) .AND. (KTYPE /=1) ) THEN
  CALL ABOR1('TPM_FFTW:EXEC_EFFTW : WRONG VALUE KTYPE')
ENDIF

IF( LD_ALL )THEN
  CALL CREATE_PLAN_FFTW(IPLAN_C2R,KTYPE,KRLEN,KFIELDS)
  IF  (JPRB == JPRD) THEN
     ZFFTP=FFTW_ALLOC_COMPLEX(INT(KCLEN/2*KFIELDS,C_SIZE_T))
  ELSE
     ZFFTP=FFTWF_ALLOC_COMPLEX(INT(KCLEN/2*KFIELDS,C_SIZE_T))
  END IF
  CALL C_F_POINTER(ZFFTP,ZFFT,[KCLEN,KFIELDS])
  IF (KTYPE==1) THEN
    DO JF=1,KFIELDS
      DO JJ=1,KCLEN
        ZFFT(JJ,JF) =PREEL(KOFF+JJ-1,JF)
      ENDDO
    ENDDO
    IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_C2R',0,ZHOOK_HANDLE2)
    IF (JPRB == JPRD) THEN
       CALL DFFTW_EXECUTE_DFT_C2R(IPLAN_C2R,ZFFT,ZFFT)
    ELSE
       CALL SFFTW_EXECUTE_DFT_C2R(IPLAN_C2R,ZFFT,ZFFT)
    END IF
    IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_C2R',1,ZHOOK_HANDLE2)
    DO JJ=1,KRLEN
      DO JF=1,KFIELDS
        PREEL(KOFF+JJ-1,JF)=ZFFT(JJ,JF)
      ENDDO
    ENDDO
  ELSE
    DO JF=1,KFIELDS
      DO JJ=1,KRLEN
        ZFFT(JJ,JF) =PREEL(KOFF+JJ-1,JF)
      ENDDO
    ENDDO
    IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_R2C',0,ZHOOK_HANDLE2)
    IF (JPRB == JPRD) THEN
       CALL DFFTW_EXECUTE_DFT_R2C(IPLAN_C2R,ZFFT,ZFFT)
    ELSE
       CALL SFFTW_EXECUTE_DFT_R2C(IPLAN_C2R,ZFFT,ZFFT)
    END IF
    IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_R2C',1,ZHOOK_HANDLE2)
    DO JJ=1,KCLEN
      DO JF=1,KFIELDS
        PREEL(KOFF+JJ-1,JF)=ZFFT(JJ,JF)/REAL(KRLEN,JPRB)
      ENDDO
    ENDDO
  ENDIF
  IF (JPRB == JPRD) THEN
     CALL FFTW_FREE(ZFFTP)
  ELSE
     CALL FFTWF_FREE(ZFFTP)
  END IF 
ELSE
  CALL CREATE_PLAN_FFTW(IPLAN_C2R1,KTYPE,KRLEN,1)
  IF (JPRB == JPRD) THEN
     ZFFT1P=FFTW_ALLOC_COMPLEX(INT(KCLEN/2,C_SIZE_T))
  ELSE
     ZFFT1P=FFTWF_ALLOC_COMPLEX(INT(KCLEN/2,C_SIZE_T))
  END IF
  CALL C_F_POINTER(ZFFT1P,ZFFT1,[KCLEN])
  IF (KTYPE==1) THEN
    DO JF=1,KFIELDS
      DO JJ=1,KCLEN
        ZFFT1(JJ) =PREEL(KOFF+JJ-1,JF)
      ENDDO
      IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_C2R',0,ZHOOK_HANDLE2)
      IF (JPRB == JPRD) THEN
         CALL DFFTW_EXECUTE_DFT_C2R(IPLAN_C2R1,ZFFT1,ZFFT1)
      ELSE
         CALL SFFTW_EXECUTE_DFT_C2R(IPLAN_C2R1,ZFFT1,ZFFT1)
      END IF
      IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_C2R',1,ZHOOK_HANDLE2)
      DO JJ=1,KRLEN
        PREEL(KOFF+JJ-1,JF)=ZFFT1(JJ)
      ENDDO
    ENDDO
  ELSE
    DO JF=1,KFIELDS
      DO JJ=1,KRLEN
        ZFFT1(JJ) =PREEL(KOFF+JJ-1,JF)
      ENDDO
      IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_R2C',0,ZHOOK_HANDLE2)
      IF (JPRB == JPRD) THEN
         CALL DFFTW_EXECUTE_DFT_R2C(IPLAN_C2R1,ZFFT1,ZFFT1)
      ELSE
         CALL SFFTW_EXECUTE_DFT_R2C(IPLAN_C2R1,ZFFT1,ZFFT1)
      END IF
      IF (LHOOK) CALL DR_HOOK('FFTW_EXECUTE_DFT_R2C',1,ZHOOK_HANDLE2)
      DO JJ=1,KCLEN
        PREEL(KOFF+JJ-1,JF)=ZFFT1(JJ)/REAL(KRLEN,JPRB)
      ENDDO
    ENDDO
  ENDIF
  IF (JPRB == JPRD) THEN
     CALL FFTW_FREE(ZFFT1P)
  ELSE
     CALL FFTWF_FREE(ZFFT1P)
  END IF
ENDIF

IF (LHOOK) CALL DR_HOOK('EXEC_EFFTW',1,ZHOOK_HANDLE)
END SUBROUTINE EXEC_EFFTW

END MODULE TPM_FFTW