!
!  $Author: tomita $
!  $Date: 2007/08/01 20:09:58 $
!  $Revision: 1.1.1.1 $
!
MODULE Fourier

!  USE ModConstants, ONLY : r8
  use typeKinds, only: r8 => Double
!                              Ifmx => MFactorFourier,  &
!                              Imxt => MTrigsFourier

  IMPLICIT NONE

  PRIVATE

  PUBLIC :: InitFFT, fft991

  INTEGER, PARAMETER :: MFactorFourier = 64 
  INTEGER :: MTrigsFourier

!  INTEGER, DIMENSION (:), ALLOCATABLE :: ifax
!  REAL (KIND=r8), DIMENSION (:), ALLOCATABLE :: trigs


CONTAINS


SUBROUTINE InitFFT(iMax, ifax, trigs)
  INTEGER,                    INTENT(IN   ) :: iMax
  INTEGER,       ALLOCATABLE, INTENT(INOUT) :: ifax(:)
  REAL(KIND=r8), ALLOCATABLE, INTENT(INOUT) :: trigs(:)
  
  MTrigsFourier  = 3*iMax/2

  ALLOCATE(ifax(MFactorFourier))
  ALLOCATE(trigs(MtrigsFourier))

  CALL fax    (iMax, 2, ifax)
  CALL fftrig (iMax, 2, trigs)

END SUBROUTINE InitFFT



!SUBROUTINE destroyFFT
!
!  IMPLICIT NONE
!
!  DEALLOCATE (ifax)
!  DEALLOCATE (trigs)
!
!END SUBROUTINE destroyFFT


SUBROUTINE fax (n, mode, ifax)

  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: n
  INTEGER, INTENT(IN   ) :: mode
  INTEGER, INTENT(INOUT) :: ifax(:)

  INTEGER :: i, l, k, ii, nn, inc, item, nfax, istop

  nn=n
  IF (ABS(mode) == 1) GO TO 10
  IF (ABS(mode) == 8) GO TO 10
  nn=n/2
  IF ((nn+nn) == n) GO TO 10
  ifax(1)=-99
  RETURN

10 CONTINUE
  k=1
20 CONTINUE
  IF (MOD(nn,4) /= 0) GO TO 30
  k=k+1
  ifax(k)=4
  nn=nn/4
  IF (nn == 1) GO TO 80
  GO TO 20
30 CONTINUE
  IF (MOD(nn,2) /= 0) GO TO 40
  k=k+1
  ifax(k)=2
  nn=nn/2
  IF (nn == 1) GO TO 80
40 CONTINUE
  IF (MOD(nn,3) /= 0) GO TO 50
  k=k+1
  ifax(k)=3
  nn=nn/3
  IF (nn == 1) GO TO 80
  GO TO 40
50 CONTINUE
  l=5
  inc=2
60 CONTINUE
  IF (MOD(nn,l) /= 0) GO TO 70
  k=k+1
  ifax(k)=l
  nn=nn/l
  IF (nn == 1) GO TO 80
  GO TO 60
70 CONTINUE
  l=l+inc
  inc=6-inc
  GO TO 60
80 CONTINUE
  ifax(1)=k-1
  nfax=ifax(1)
  IF (nfax /= 1) THEN
     DO ii=2,nfax
        istop=nfax+2-ii
        DO i=2,istop
           IF (ifax(i+1) >= ifax(i)) CYCLE
           item=ifax(i)
           ifax(i)=ifax(i+1)
           ifax(i+1)=item
        END DO
     END DO
  END IF

END SUBROUTINE fax


SUBROUTINE fftrig (n, mode, trigs)

  IMPLICIT NONE

  INTEGER,       INTENT(IN   ) :: n
  INTEGER,       INTENT(IN   ) :: mode
  REAL(KIND=r8), INTENT(INOUT) :: trigs(:)

  INTEGER :: i, l, la, nh, nn, imode

  REAL (KIND=r8) :: pi, del, angle

  pi=2.0_r8*ASIN(1.0_r8)
  imode=ABS(mode)
  nn=n
  IF (imode > 1 .AND. imode < 6) nn=n/2
  del=(pi+pi)/REAL(nn,r8)
  l=nn+nn
  DO i=1,l,2
     angle=0.5_r8*REAL(i-1,r8)*del
     trigs(i)=COS(angle)
     trigs(i+1)=SIN(angle)
  END DO
  IF (imode == 1) RETURN
  IF (imode == 8) RETURN
  del=0.5_r8*del
  nh=(nn+1)/2
  l=nh+nh
  la=nn+nn
  DO i=1,l,2
     angle=0.5_r8*REAL(i-1,r8)*del
     trigs(la+i)=COS(angle)
     trigs(la+i+1)=SIN(angle)
  END DO
  IF (imode <= 3) RETURN
  del=0.5_r8*del
  la=la+nn
  IF (mode /= 5) THEN
     DO i=2,nn
        angle=REAL(i-1,r8)*del
        trigs(la+i)=2.0_r8*SIN(angle)
     END DO
  ELSE
     del=0.5_r8*del
     DO i=2,n
        angle=REAL(i-1,r8)*del
        trigs(la+i)=SIN(angle)
     END DO
  END IF

END SUBROUTINE fftrig


SUBROUTINE fft991 (a, work, inc, jump, n, lot, ifax, trigs, isign)

  ! fft991 :transforms grid point representations of fields into their
  !         fourier representations at all model levels and vice-versa.
  !
  !        i=sqrt(-1)
  !        0<x<2*pai
  !        max=n/2
  !     f(x)=f(0)+f(  1)*exp(-i*  1*x)+f(-  1)*exp(+i*  1*x)
  !              +f(  2)*exp(-i*  2*x)+f(-  2)*exp(+i*  2*x)
  !              +........................
  !              +f(  m)*exp(-i*  m*x)+f(-  m)*exp(+i*  m*x)
  !              +........................
  !              +f(max)*exp(-i*max*x)+f(-max)*exp(+i*max*x)
  !-----------------------------------------------------------------------
  !..  a    (jump,lot)   isign=+1
  !                 input array  for m=0,1,2,....,n/2
  !                      a(2*m+1)=real(f(m))
  !                      a(2*m+2)=imag(f(m))
  !                      a(  1)=cos(m= 0)   ,a(2)=0.
  !                      a(  3)=cos(m= 1)   ,a(4)=sin(m=1)
  !                      a(  5)=cos(m= 2)   ,a(6)=sin(m=2)
  !                      .................................................
  !                      a(n-1)=cos(m=n/2-1),a(n)=sin(m=n/2-1)
  !                      a(n+1)=cos(m=n/2)
  !                 output array
  !                      grid point values
  !                      a(n+1)=0.,......a(jump)=0.
  !..  a    (jump,lot)  isign=-1
  !                 input array
  !                      grid point values
  !                      a(n+1)=0.,......a(jump)=0.
  !                 output array  for m=0,1,2,....,n/2
  !                      a(2*m+1)=real(f(m))
  !                      a(2*m+2)=imag(f(m))
  !                      a(  1)=cos(m= 0)   ,a(2)=0.
  !                      a(  3)=cos(m= 1)   ,a(4)=sin(m= 1)
  !                      a(  5)=cos(m= 2)   ,a(6)=sin(m= 2)
  !                      .................................................
  !                      a(n-1)=cos(m=n/2-1),a(n)=sin(m=n/2-1)
  !                      a(n+1)=cos(m=n/2)
  !..  work (jump,lot)     work  array
  !..  trigs(.gt.1.5*Imax) sin,cos coefficient
  !..  ifax (*)            factors of Imax
  !..  inc                 increment of data spacing
  !..  jump                jump.ge.n+1
  !..  n                   number of grids to be trasformed
  !..  lot                 number of lines transformed simultaneously
  !..  isign=+1            wave to grid transform
  !          -1            grid to wave transform
  !-----------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER, INTENT(IN) :: inc, jump, n, lot, isign

  INTEGER,        DIMENSION(:), INTENT(IN   ) :: ifax
  REAL (KIND=r8), DIMENSION(:), INTENT(IN   ) :: trigs
  REAL (KIND=r8), DIMENSION(:), INTENT(INOUT) :: a
  REAL (KIND=r8), DIMENSION(:), INTENT(INOUT) :: work

  INTEGER :: i, j, k, l, m, ia, ib, la, nx, nh
  INTEGER :: ink, igo, nfax, ibase, jbase

  nfax=ifax(1)
  nx=n+1
  nh=n/2
  ink=inc+inc
  IF (isign == 1) THEN
     CALL fft99a (a, work, trigs, inc, jump, n, lot)
     igo=60
  ELSE
     igo=50
     IF (MOD(nfax,2) /= 1)THEN
        ibase=1
        jbase=1
        DO l=1,lot
           i=ibase
           j=jbase
!cdir nodep
           DO m=1,n
              work(j)=a(i)
              i=i+inc
              j=j+1
           END DO
           ibase=ibase+jump
           jbase=jbase+nx
        END DO
        igo=60
     END IF
  END IF
  ia=1
  la=1
  DO k=1,nfax
     IF (igo == 60) THEN
        CALL vpassm (work(1:), work(2:), a(ia:), a(ia+inc:), trigs, &
                     2, ink, nx, jump, lot, nh, ifax(k+1), la)
        igo=50
     ELSE
        CALL vpassm (a(ia:), a(ia+inc:), work(1:), work(2:), trigs,&
                     ink, 2, jump, nx, lot, nh, ifax(k+1), la)
        igo=60
     END IF
     la=la*ifax(k+1)
  END DO
  IF (isign == -1) THEN
     CALL fft99b (work, a, trigs, inc, jump, n, lot)
  ELSE
     IF (MOD(nfax,2) /= 1)THEN
        ibase=1
        jbase=1
        DO l=1,lot
           i=ibase
           j=jbase
!cdir nodep
           DO m=1,n
              a(j)=work(i)
              i=i+1
              j=j+inc
           END DO
           ibase=ibase+nx
           jbase=jbase+jump
        END DO
     END IF
     ib=n*inc+1
!cdir nodep
     DO l=1,lot
        a(ib)=0.0_r8
        a(ib+inc)=0.0_r8
        ib=ib+jump
     END DO
  END IF

END SUBROUTINE fft991


SUBROUTINE fft99a (a, work, trigs, inc, jump, n, lot)

  IMPLICIT NONE

  INTEGER, INTENT(IN) :: inc, jump, n, lot
  REAL (KIND=r8), DIMENSION(:), INTENT(IN   ) :: a
  REAL (KIND=r8), DIMENSION(:), INTENT(IN   ) :: trigs
  REAL (KIND=r8), DIMENSION(:), INTENT(  OUT) :: work

  INTEGER :: l, k, ia, ib, ja, jb, nh, nx, ink
  INTEGER :: iabase, ibbase, jabase, jbbase

  REAL (KIND=r8) :: c, s

  nh=n/2
  nx=n+1
  ink=inc+inc
  ia=1
  ib=n*inc+1
  ja=1
  jb=2
!cdir nodep
  DO l=1,lot
     work(ja)=a(ia)+a(ib)
     work(jb)=a(ia)-a(ib)
     ia=ia+jump
     ib=ib+jump
     ja=ja+nx
     jb=jb+nx
  END DO
  iabase=2*inc+1
  ibbase=(n-2)*inc+1
  jabase=3
  jbbase=n-1
  DO k=3,nh,2
     ia=iabase
     ib=ibbase
     ja=jabase
     jb=jbbase
     c=trigs(n+k)
     s=trigs(n+k+1)
!cdir nodep
     DO l=1,lot
        work(ja)=(a(ia)+a(ib))- &
                 (s*(a(ia)-a(ib))+c*(a(ia+inc)+a(ib+inc)))
        work(jb)=(a(ia)+a(ib))+ &
                 (s*(a(ia)-a(ib))+c*(a(ia+inc)+a(ib+inc)))
        work(ja+1)=(c*(a(ia)-a(ib))-s*(a(ia+inc)+a(ib+inc)))+ &
                      (a(ia+inc)-a(ib+inc))
        work(jb+1)=(c*(a(ia)-a(ib))-s*(a(ia+inc)+a(ib+inc)))- &
                      (a(ia+inc)-a(ib+inc))
        ia=ia+jump
        ib=ib+jump
        ja=ja+nx
        jb=jb+nx
     END DO
     iabase=iabase+ink
     ibbase=ibbase-ink
     jabase=jabase+2
     jbbase=jbbase-2
  END DO
  IF (iabase == ibbase) THEN
     ia=iabase
     ja=jabase
!cdir nodep
     DO l=1,lot
        work(ja)=2.0_r8*a(ia)
        work(ja+1)=-2.0_r8*a(ia+inc)
        ia=ia+jump
        ja=ja+nx
     END DO
  END IF

END SUBROUTINE fft99a


SUBROUTINE fft99b (work, a, trigs, inc, jump, n, lot)

  IMPLICIT NONE

  INTEGER, INTENT(IN) :: inc, jump, n, lot

  REAL (KIND=r8), DIMENSION(:), INTENT(IN   ) :: work
  REAL (KIND=r8), DIMENSION(:), INTENT(IN   ) :: trigs
  REAL (KIND=r8), DIMENSION(:), INTENT(OUT  ) :: a

  INTEGER :: l, k, ia, ib, ja, jb, nh, nx, ink
  INTEGER :: iabase, ibbase, jabase, jbbase

  REAL (KIND=r8) :: c, s, scale

  nh=n/2
  nx=n+1
  ink=inc+inc
  scale=1.0_r8/REAL(n,r8)
  ia=1
  ib=2
  ja=1
  jb=n*inc+1
!cdir nodep
  DO l=1,lot
     a(ja)=scale*(work(ia)+work(ib))
     a(jb)=scale*(work(ia)-work(ib))
     a(ja+inc)=0.0_r8
     a(jb+inc)=0.0_r8
     ia=ia+nx
     ib=ib+nx
     ja=ja+jump
     jb=jb+jump
  END DO
  scale=0.5_r8*scale
  iabase=3
  ibbase=n-1
  jabase=2*inc+1
  jbbase=(n-2)*inc+1
  DO k=3,nh,2
     ia=iabase
     ib=ibbase
     ja=jabase
     jb=jbbase
     c=trigs(n+k)
     s=trigs(n+k+1)
!cdir nodep
     DO l=1,lot
        a(ja)=scale*((work(ia)+work(ib)) &
             +(c*(work(ia+1)+work(ib+1))+s*(work(ia)-work(ib))))
        a(jb)=scale*((work(ia)+work(ib)) &
             -(c*(work(ia+1)+work(ib+1))+s*(work(ia)-work(ib))))
        a(ja+inc)=scale*((c*(work(ia)-work(ib))-s*(work(ia+1)+ &
             work(ib+1)))+(work(ib+1)-work(ia+1)))
        a(jb+inc)=scale*((c*(work(ia)-work(ib))-s*(work(ia+1)+ &
             work(ib+1)))-(work(ib+1)-work(ia+1)))
        ia=ia+nx
        ib=ib+nx
        ja=ja+jump
        jb=jb+jump
     END DO
     iabase=iabase+2
     ibbase=ibbase-2
     jabase=jabase+ink
     jbbase=jbbase-ink
  END DO
  IF (iabase == ibbase)THEN
     ia=iabase
     ja=jabase
     scale=2.0_r8*scale
!cdir nodep
     DO l=1,lot
        a(ja)=scale*work(ia)
        a(ja+inc)=-scale*work(ia+1)
        ia=ia+nx
        ja=ja+jump
     END DO
  END IF

END SUBROUTINE fft99b


SUBROUTINE vpassm (a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la)

  IMPLICIT NONE

  INTEGER, INTENT(IN) :: inc1, inc2, inc3, inc4, lot, n, ifac, la
  
  REAL (KIND=r8), DIMENSION (:), INTENT(IN   )  :: a, b
  REAL (KIND=r8), DIMENSION (:), INTENT(OUT  ) :: c, d
  REAL (KIND=r8), DIMENSION (:), INTENT(IN   ) :: trigs

  INTEGER :: i, j, k, l, m
  INTEGER :: ia, ja, ib, jb, kb, ic, jc, kc, id, jd, kd, ie, je, ke
  INTEGER :: igo, ijk, la1
  INTEGER :: iink, jink, jump, ibase, jbase

  REAL (KIND=r8) :: c1, s1,c2, s2, c3, s3, c4, s4
  REAL (KIND=r8) :: wka, wkb, wksina, wksinb, wkaacp, wkbacp, wkaacm, wkbacm

  INTEGER, SAVE :: ifr=0

  REAL (KIND=r8), SAVE :: radi,sin60,sin36,sin72,cos36,cos72

  IF (ifr == 0) THEN
     radi=ATAN(1.0_r8)/45.0_r8
     sin60=SIN(60.0_r8*radi)
     sin36=SIN(36.0_r8*radi)
     sin72=SIN(72.0_r8*radi)
     cos36=COS(36.0_r8*radi)
     cos72=COS(72.0_r8*radi)
     ifr=1
  END IF
  m=n/ifac
  iink=m*inc1
  jink=la*inc2
  jump=(ifac-1)*jink
  ibase=0
  jbase=0
  igo=ifac-1
  IF (igo == 1) THEN
     ia=1
     ja=1
     ib=ia+iink
     jb=ja+jink
     DO l=1,la
        i=ibase
        j=jbase
!cdir nodep
        DO ijk=1,lot
           c(ja+j)=a(ia+i)+a(ib+i)
           d(ja+j)=b(ia+i)+b(ib+i)
           c(jb+j)=a(ia+i)-a(ib+i)
           d(jb+j)=b(ia+i)-b(ib+i)
           i=i+inc3
           j=j+inc4
        END DO
        ibase=ibase+inc1
        jbase=jbase+inc2
     END DO
     IF (la == m) RETURN
     la1=la+1
     jbase=jbase+jump
     DO k=la1,m,la
        kb=k+k-2
        c1=trigs(kb+1)
        s1=trigs(kb+2)
        DO l=1,la
           i=ibase
           j=jbase
!cdir nodep
           DO ijk=1,lot
              wka=a(ia+i)-a(ib+i)
              wkb=b(ia+i)-b(ib+i)
              c(ja+j)=a(ia+i)+a(ib+i)
              d(ja+j)=b(ia+i)+b(ib+i)
              c(jb+j)=c1*wka-s1*wkb
              d(jb+j)=s1*wka+c1*wkb
              i=i+inc3
              j=j+inc4
           END DO
           ibase=ibase+inc1
           jbase=jbase+inc2
        END DO
        jbase=jbase+jump
     END DO
  ELSEIF (igo == 2) THEN
     ia=1
     ja=1
     ib=ia+iink
     jb=ja+jink
     ic=ib+iink
     jc=jb+jink
     DO l=1,la
        i=ibase
        j=jbase
!cdir  nodep
        DO ijk=1,lot
           wka=a(ib+i)+a(ic+i)
           wkb=b(ib+i)+b(ic+i)
           wksina=sin60*(a(ib+i)-a(ic+i))
           wksinb=sin60*(b(ib+i)-b(ic+i))
           c(ja+j)=a(ia+i)+wka
           d(ja+j)=b(ia+i)+wkb
           c(jb+j)=(a(ia+i)-0.5_r8*wka)-wksinb
           c(jc+j)=(a(ia+i)-0.5_r8*wka)+wksinb
           d(jb+j)=(b(ia+i)-0.5_r8*wkb)+wksina
           d(jc+j)=(b(ia+i)-0.5_r8*wkb)-wksina
           i=i+inc3
           j=j+inc4
        END DO
        ibase=ibase+inc1
        jbase=jbase+inc2
     END DO
     IF (la == m) RETURN
     la1=la+1
     jbase=jbase+jump
     DO k=la1,m,la
        kb=k+k-2
        kc=kb+kb
        c1=trigs(kb+1)
        s1=trigs(kb+2)
        c2=trigs(kc+1)
        s2=trigs(kc+2)
        DO l=1,la
           i=ibase
           j=jbase
!cdir nodep
           DO ijk=1,lot
              wka=a(ib+i)+a(ic+i)
              wkb=b(ib+i)+b(ic+i)
              wksina=sin60*(a(ib+i)-a(ic+i))
              wksinb=sin60*(b(ib+i)-b(ic+i))
              c(ja+j)=a(ia+i)+wka
              d(ja+j)=b(ia+i)+wkb
              c(jb+j)=c1*((a(ia+i)-0.5_r8*wka)-wksinb) &
                     -s1*((b(ia+i)-0.5_r8*wkb)+wksina)
              d(jb+j)=s1*((a(ia+i)-0.5_r8*wka)-wksinb) &
                     +c1*((b(ia+i)-0.5_r8*wkb)+wksina)
              c(jc+j)=c2*((a(ia+i)-0.5_r8*wka)+wksinb) &
                     -s2*((b(ia+i)-0.5_r8*wkb)-wksina)
              d(jc+j)=s2*((a(ia+i)-0.5_r8*wka)+wksinb) &
                     +c2*((b(ia+i)-0.5_r8*wkb)-wksina)
              i=i+inc3
              j=j+inc4
           END DO
           ibase=ibase+inc1
           jbase=jbase+inc2
        END DO
        jbase=jbase+jump
     END DO
  ELSEIF (igo == 3) THEN
     ia=1
     ja=1
     ib=ia+iink
     jb=ja+jink
     ic=ib+iink
     jc=jb+jink
     id=ic+iink
     jd=jc+jink
     DO l=1,la
        i=ibase
        j=jbase
!cdir nodep
        DO ijk=1,lot
           wkaacp=a(ia+i)+a(ic+i)
           wkbacp=b(ia+i)+b(ic+i)
           wkaacm=a(ia+i)-a(ic+i)
           wkbacm=b(ia+i)-b(ic+i)
           c(ja+j)=wkaacp+(a(ib+i)+a(id+i))
           c(jc+j)=wkaacp-(a(ib+i)+a(id+i))
           d(ja+j)=wkbacp+(b(ib+i)+b(id+i))
           d(jc+j)=wkbacp-(b(ib+i)+b(id+i))
           c(jb+j)=wkaacm-(b(ib+i)-b(id+i))
           c(jd+j)=wkaacm+(b(ib+i)-b(id+i))
           d(jb+j)=wkbacm+(a(ib+i)-a(id+i))
           d(jd+j)=wkbacm-(a(ib+i)-a(id+i))
           i=i+inc3
           j=j+inc4
        END DO
        ibase=ibase+inc1
        jbase=jbase+inc2
     END DO
     IF (la == m) RETURN
     la1=la+1
     jbase=jbase+jump
     DO k=la1,m,la
        kb=k+k-2
        kc=kb+kb
        kd=kc+kb
        c1=trigs(kb+1)
        s1=trigs(kb+2)
        c2=trigs(kc+1)
        s2=trigs(kc+2)
        c3=trigs(kd+1)
        s3=trigs(kd+2)
        DO l=1,la
           i=ibase
           j=jbase
!cdir nodep
           DO ijk=1,lot
              wkaacp=a(ia+i)+a(ic+i)
              wkbacp=b(ia+i)+b(ic+i)
              wkaacm=a(ia+i)-a(ic+i)
              wkbacm=b(ia+i)-b(ic+i)
              c(ja+j)=wkaacp+(a(ib+i)+a(id+i))
              d(ja+j)=wkbacp+(b(ib+i)+b(id+i))
              c(jc+j)= c2*(wkaacp-(a(ib+i)+a(id+i))) &
                      -s2*(wkbacp-(b(ib+i)+b(id+i))) 
              d(jc+j)= s2*(wkaacp-(a(ib+i)+a(id+i))) &
                      +c2*(wkbacp-(b(ib+i)+b(id+i)))
              c(jb+j)= c1*(wkaacm-(b(ib+i)-b(id+i))) &
                      -s1*(wkbacm+(a(ib+i)-a(id+i)))
              d(jb+j)= s1*(wkaacm-(b(ib+i)-b(id+i))) &
                      +c1*(wkbacm+(a(ib+i)-a(id+i)))
              c(jd+j)= c3*(wkaacm+(b(ib+i)-b(id+i))) &
                      -s3*(wkbacm-(a(ib+i)-a(id+i)))
              d(jd+j)= s3*(wkaacm+(b(ib+i)-b(id+i))) &
                      +c3*(wkbacm-(a(ib+i)-a(id+i)))
              i=i+inc3
              j=j+inc4
           END DO
           ibase=ibase+inc1
           jbase=jbase+inc2
        END DO
        jbase=jbase+jump
     END DO
  ELSEIF (igo == 4) THEN
     ia=1
     ja=1
     ib=ia+iink
     jb=ja+jink
     ic=ib+iink
     jc=jb+jink
     id=ic+iink
     jd=jc+jink
     ie=id+iink
     je=jd+jink
     DO l=1,la
        i=ibase
        j=jbase
!cdir nodep
        DO ijk=1,lot
           c(ja+j)=a(ia+i)+(a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i))
           d(ja+j)=b(ia+i)+(b(ib+i)+b(ie+i))+(b(ic+i)+b(id+i))
           c(jb+j)=(a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+ &
                a(id+i)))-(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i)))
           c(je+j)=(a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+ &
                a(id+i)))+(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i)))
           d(jb+j)=(b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+ &
                b(id+i)))+(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i)))
           d(je+j)=(b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+ &
                b(id+i)))-(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i)))
           c(jc+j)=(a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+ &
                a(id+i)))-(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i)))
           c(jd+j)=(a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+ &
                a(id+i)))+(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i)))
           d(jc+j)=(b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+ &
                b(id+i)))+(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i)))
           d(jd+j)=(b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+ &
                b(id+i)))-(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i)))
           i=i+inc3
           j=j+inc4
        END DO
        ibase=ibase+inc1
        jbase=jbase+inc2
     END DO
     IF (la == m) RETURN
     la1=la+1
     jbase=jbase+jump
     DO k=la1,m,la
        kb=k+k-2
        kc=kb+kb
        kd=kc+kb
        ke=kd+kb
        c1=trigs(kb+1)
        s1=trigs(kb+2)
        c2=trigs(kc+1)
        s2=trigs(kc+2)
        c3=trigs(kd+1)
        s3=trigs(kd+2)
        c4=trigs(ke+1)
        s4=trigs(ke+2)
        DO l=1,la
           i=ibase
           j=jbase
!cdir nodep
           DO ijk=1,lot
              c(ja+j)=a(ia+i)+(a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i))
              d(ja+j)=b(ia+i)+(b(ib+i)+b(ie+i))+(b(ic+i)+b(id+i))
              c(jb+j)=c1*((a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36* &
                   (a(ic+i)+a(id+i)))-(sin72*(b(ib+i)-b(ie+i)) &
                   +sin36*(b(ic+i)-b(id+i)))) &
                   -s1*((b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36* &
                   (b(ic+i)+b(id+i)))+(sin72*(a(ib+i)-a(ie+i)) &
                   +sin36*(a(ic+i)-a(id+i))))
              d(jb+j)=s1*((a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36* &
                   (a(ic+i)+a(id+i)))-(sin72*(b(ib+i)-b(ie+i)) &
                   +sin36*(b(ic+i)-b(id+i)))) &
                   +c1*((b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36* &
                   (b(ic+i)+b(id+i)))+(sin72*(a(ib+i)-a(ie+i)) &
                   +sin36*(a(ic+i)-a(id+i))))
              c(je+j)=c4*((a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36* &
                   (a(ic+i)+a(id+i)))+(sin72*(b(ib+i)-b(ie+i)) &
                   +sin36*(b(ic+i)-b(id+i)))) &
                   -s4*((b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36* &
                   (b(ic+i)+b(id+i)))-(sin72*(a(ib+i)-a(ie+i))+ &
                   sin36*(a(ic+i)-a(id+i))))
              d(je+j)=s4*((a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36* &
                   (a(ic+i)+a(id+i)))+(sin72*(b(ib+i)-b(ie+i))+ &
                   sin36*(b(ic+i)-b(id+i))))+c4*((b(ia+i)+ &
                   cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i))) &
                   -(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i))))
              c(jc+j)=c2*((a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72* &
                   (a(ic+i)+a(id+i)))-(sin36*(b(ib+i)-b(ie+i))- &
                   sin72*(b(ic+i)-b(id+i))))-s2*((b(ia+i)-cos36* &
                   (b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i)))+(sin36 &
                   *(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i)))) 
              d(jc+j)=s2*((a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72* &
                   (a(ic+i)+a(id+i)))-(sin36*(b(ib+i)-b(ie+i))- &
                   sin72*(b(ic+i)-b(id+i))))+c2*((b(ia+i)-cos36* &
                   (b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i)))+(sin36 &
                   *(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i))))
              c(jd+j)=c3*((a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72* &
                   (a(ic+i)+a(id+i)))+(sin36*(b(ib+i)-b(ie+i))- &
                   sin72*(b(ic+i)-b(id+i))))-s3*((b(ia+i)-cos36* &
                   (b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i)))-(sin36 &
                   *(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i))))
              d(jd+j)=s3*((a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72* &
                   (a(ic+i)+a(id+i)))+(sin36*(b(ib+i)-b(ie+i))-sin72 &
                   *(b(ic+i)-b(id+i))))+c3*((b(ia+i)-cos36*(b(ib+i)+ &
                   b(ie+i))+cos72*(b(ic+i)+b(id+i)))-(sin36* &
                   (a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i))))
              i=i+inc3
              j=j+inc4
           END DO
           ibase=ibase+inc1
           jbase=jbase+inc2
        END DO
        jbase=jbase+jump
     END DO
  ENDIF

END SUBROUTINE vpassm


END MODULE Fourier