MODULE m_Expansion

USE m_Parameters, ONLY: ki8, kr8, Imax, Jmax, Imax2, Jmaxh, Jmaxh1

USE m_Arrays, ONLY: guns, guna, gvns, gvna, gpns, gpna, &
                    uh, vh, ph, w, un, vn, pn, ur, vr, pr, &
                    wunn, wuns, wvnn, wvns, wpnn, wpns

USE m_Fourier, ONLY: fft991

IMPLICIT NONE

PRIVATE

PUBLIC :: Coefficients_Fourier, Coefficients_Hough, &
          Recomposition_Fourier, Recomposition_Hough


CONTAINS


SUBROUTINE Coefficients_Fourier
 
  IMPLICIT NONE
 
  ! Perfoms Fourier Transform in Complex Domain.

  INTEGER (KIND=ki8) :: j, jj

  ! Fourier Coefficients by Symmetry

  DO j=1,Jmaxh
    jj=Jmax-j+1
    guns(1:Imax,j)=un(:,j)+un(:,jj)
    guna(1:Imax,j)=un(:,j)-un(:,jj)
    gvns(1:Imax,j)=vn(:,j)+vn(:,jj)
    gvna(1:Imax,j)=vn(:,j)-vn(:,jj)
    gpns(1:Imax,j)=pn(:,j)+pn(:,jj)
    gpna(1:Imax,j)=pn(:,j)-pn(:,jj)
    guns(Imax+1:Imax2,j)=0.0
    guna(Imax+1:Imax2,j)=0.0
    gvns(Imax+1:Imax2,j)=0.0
    gvna(Imax+1:Imax2,j)=0.0
    gpns(Imax+1:Imax2,j)=0.0
    gpna(Imax+1:Imax2,j)=0.0
  END DO

  CALL Grid_to_Fourier (guns)
  CALL Grid_to_Fourier (gvns)
  CALL Grid_to_Fourier (gpns)
  CALL Grid_to_Fourier (guna)
  CALL Grid_to_Fourier (gvna)
  CALL Grid_to_Fourier (gpna)

END SUBROUTINE Coefficients_Fourier


SUBROUTINE Coefficients_Hough (k, l, isi, cft)
 
  IMPLICIT NONE
 
  ! Performs Hough Function Decomposition

  INTEGER (KIND=ki8), INTENT (IN) :: k, l, isi

  COMPLEX (KIND=kr8), INTENT (OUT) :: cft

  INTEGER (KIND=ki8) :: ks, kc, j

  REAL (KIND=kr8) :: cfr, cfi

  REAL (KIND=kr8), SAVE :: epsm=2.0**(-39)

  cft=CMPLX(0.0,0.0)
  ks=2*k
  kc=ks-1
  cfr=0.0
  cfi=0.0
  IF (isi == 1) THEN
    DO j=1,Jmaxh
      cfr=cfr+(guns(kc,j)*uh(k,l,isi,j)-gvna(ks,j)*vh(k,l,isi,j)+gpns(kc,j)*ph(k,l,isi,j))*w(j)
      cfi=cfi+(guns(ks,j)*uh(k,l,isi,j)+gvna(kc,j)*vh(k,l,isi,j)+gpns(ks,j)*ph(k,l,isi,j))*w(j)
    END DO
  ELSE
    DO j=1,Jmaxh
      cfr=cfr+(guna(kc,j)*uh(k,l,isi,j)-gvns(ks,j)*vh(k,l,isi,j)+gpna(kc,j)*ph(k,l,isi,j))*w(j)
      cfi=cfi+(guna(ks,j)*uh(k,l,isi,j)+gvns(kc,j)*vh(k,l,isi,j)+gpna(ks,j)*ph(k,l,isi,j))*w(j)
    END DO
  END IF
  cft=CMPLX(cfr,cfi)
  IF (ABS(cft) <= epsm) cft=CMPLX(0.0,0.0)

END SUBROUTINE Coefficients_Hough


SUBROUTINE Recomposition_Fourier (nc)
 
  IMPLICIT NONE
 
  ! Recompose the Carrying Variables over Fourier and Symmetry.

  INTEGER (KIND=ki8), INTENT (IN) :: nc

  INTEGER (KIND=ki8) :: j, ji, jj

  CALL Fourier_to_Grid (wunn, ur(:,1:Jmaxh,nc))
  CALL Fourier_to_Grid (wuns, ur(:,Jmaxh1:Jmax,nc))
  CALL Fourier_to_Grid (wvnn, vr(:,1:Jmaxh,nc))
  CALL Fourier_to_Grid (wvns, vr(:,Jmaxh1:Jmax,nc))
  CALL Fourier_to_Grid (wpnn, pr(:,1:Jmaxh,nc))
  CALL Fourier_to_Grid (wpns, pr(:,Jmaxh1:Jmax,nc))

  DO j=1,Jmaxh
    ji=Jmaxh-j+1
    jj=Jmaxh+j
    wunn(1:Imax,j) =ur(:,j,nc)+ur(:,jj,nc)
    wuns(1:Imax,ji)=ur(:,j,nc)-ur(:,jj,nc)
    wvnn(1:Imax,j) =vr(:,j,nc)+vr(:,jj,nc)
    wvns(1:Imax,ji)=vr(:,j,nc)-vr(:,jj,nc)
    wpnn(1:Imax,j) =pr(:,j,nc)+pr(:,jj,nc)
    wpns(1:Imax,ji)=pr(:,j,nc)-pr(:,jj,nc)
  END DO

  DO j=1,Jmaxh
    jj=Jmaxh+j
    ur(:,j,nc) =wunn(1:Imax,j)
    ur(:,jj,nc)=wuns(1:Imax,j)
    vr(:,j,nc) =wvnn(1:Imax,j)
    vr(:,jj,nc)=wvns(1:Imax,j)
    pr(:,j,nc) =wpnn(1:Imax,j)
    pr(:,jj,nc)=wpns(1:Imax,j)
  END DO

END SUBROUTINE Recomposition_Fourier


SUBROUTINE Recomposition_Hough (k, l, isi, cft)
 
  IMPLICIT NONE

  ! Recomposes the Carrying Variables Over Hough Functions.

  INTEGER (KIND=ki8), INTENT (IN) :: k, l, isi

  COMPLEX (KIND=kr8), INTENT (IN) :: cft

  INTEGER (KIND=ki8) :: ks, kc

  REAL (KIND=kr8) :: cftr, cfti

  ks=2*k
  kc=ks-1
  cftr=REAL(cft)
  cfti=AIMAG(cft)
  IF (isi == 1) THEN
    wunn(kc,:)=wunn(kc,:)+cftr*uh(k,l,isi,:)
    wunn(ks,:)=wunn(ks,:)+cfti*uh(k,l,isi,:)
    wvns(kc,:)=wvns(kc,:)+cfti*vh(k,l,isi,:)
    wvns(ks,:)=wvns(ks,:)-cftr*vh(k,l,isi,:)
    wpnn(kc,:)=wpnn(kc,:)+cftr*ph(k,l,isi,:)
    wpnn(ks,:)=wpnn(ks,:)+cfti*ph(k,l,isi,:)
  ELSE
    wuns(kc,:)=wuns(kc,:)+cftr*uh(k,l,isi,:)
    wuns(ks,:)=wuns(ks,:)+cfti*uh(k,l,isi,:)
    wvnn(kc,:)=wvnn(kc,:)+cfti*vh(k,l,isi,:)
    wvnn(ks,:)=wvnn(ks,:)-cftr*vh(k,l,isi,:)
    wpns(kc,:)=wpns(kc,:)+cftr*ph(k,l,isi,:)
    wpns(ks,:)=wpns(ks,:)+cfti*ph(k,l,isi,:)
  END IF

END SUBROUTINE Recomposition_Hough


SUBROUTINE Grid_to_Fourier (fin)
 
  IMPLICIT NONE

  ! General Grid to Spectral Fast Fourier Transform

  REAL (KIND=kr8), DIMENSION (Imax2,Jmaxh), INTENT (IN OUT) :: fin

  REAL (KIND=kr8), DIMENSION (Imax2) :: work

  INTEGER (KIND=ki8) :: inc, lot, isgn, j

  inc=1
  lot=1
  isgn=-1

  DO j=1,Jmaxh
    work=0.0
    CALL fft991 (fin(:,j), work, inc, Imax2, Imax, lot, isgn)
  END DO

END SUBROUTINE Grid_to_Fourier


SUBROUTINE Fourier_to_Grid (fin, fout)
 
  IMPLICIT NONE
 
  ! General Spectral to Grid transform

  REAL (KIND=kr8), DIMENSION (Imax2,Jmaxh), INTENT (IN OUT) :: fin

  REAL (KIND=kr8), DIMENSION (Imax,Jmaxh), INTENT (OUT) :: fout

  REAL (KIND=kr8), DIMENSION (Imax2) :: work

  INTEGER (KIND=ki8) :: inc, lot, isgn, j

  inc=1
  lot=1
  isgn=1

  DO j=1,Jmaxh
    work=0.0
    CALL fft991 (fin(:,j), work, inc, Imax2, Imax, lot, isgn)
  END DO

  fout=fin(1:Imax,:)

END SUBROUTINE Fourier_to_Grid


END MODULE m_Expansion