subroutine kparfft(f,iabs) ! ! Changing code to default complex notation. bdorland ! ! calculates i * kpar * f if iabs=-1 ! calculates kpar * f if iabs=0 ! calculates |kpar|* f if iabs=1 ! calculates -i * |kpar|* f if iabs=2 ! and f/(1+s_par**4*kpar**4) if iabs=13 ! use itg_data implicit none integer i, j, k, l, m, n, iabs, imode complex, dimension(:,:,:,:) :: f ! complex, dimension(:,:), pointer, save :: a, b ! complex, dimension(:), pointer, save :: c, d complex, dimension(ld*nconz,md*nd*nspecies) :: a, b complex, dimension(ld*nconz) :: c, d real, dimension(:), allocatable :: work, ttablekp real scale real sgpar integer isign,nwork,ntablekp interface subroutine xmcfft(isign, n, m, scale, x, inc1x, inc2x, & y, inc1y, inc2y, table, ntable, work, nwork) integer, intent (in) :: isign, n, m, inc1x, inc2x, inc1y, inc2y complex, dimension(inc2x, m) :: x complex, dimension(inc2y, m) :: y real, dimension(ntable) :: table real, dimension(nwork) :: work end subroutine xmcfft end interface nwork=4*ld*md*nd*nspecies allocate (work(nwork)) ! Is this necessary? if(iabs <= 0) then sgpar=gpar else sgpar=abs(gpar) endif ! ! easiest to use switch for iperiod=3 ! ! if(.not.associated(a)) allocate (a(ld*nconz,md*nd*nspecies), b(ld*nconz,md*nd*nspecies)) if (iperiod == 3) then ! ! do ndom FFTs, of length scon(ndom), with ncon(ndom) modes ! ! if(.not.associated(c)) allocate (c(ld*nconz), d(ld*nconz)) allocate (ttablekp(101+2*ld*(nconz+1))) do j=1,ndom ntablekp=101+2*ld*(nconz+1) scale=1.0/sqrt(float(ldb*scon(j))) do k=1,ntablekp ttablekp(k)=tablekp(j,k) enddo do k=1,ncon(j)*scon(j) do i=1,nspecies m=mmcon(j,k) n=nncon(j,k) imode=mcon(m,n)+(i-1)*ncon(j) do l=1,ldb a(l+ldb*(conpos(m,n)-1),imode)=f(l,m,n,i) enddo enddo enddo isign=-1 ! ! if number of transforms too small, use cfft ! if (ncon(j)*nspecies < 2) then do l=1,ldb*scon(j) c(l)=a(l,1) enddo call cfft(isign,ldb*scon(j),scale,c,1,d,1,ttablekp,ntablekp,work,nwork) do l=1,ldb*scon(j) b(l,1)=d(l) enddo else call xmcfft(isign,ldb*scon(j),ncon(j)*nspecies,scale,a,1,& size(a, 1),b,1,size(b, 1),ttablekp,ntablekp,work,nwork) endif if (iabs == -1) then do imode=1,ncon(j)*nspecies ! for some reason, CRAY tries to vectorize n loop by default: !cfpp$ select (vector) do l=1,ldb*scon(j) a(l,imode)=ii*kpar(l,j)*b(l,imode) enddo enddo else if (iabs == 0) then do imode=1,ncon(j)*nspecies !cfpp$ select (vector) do l=1,ldb*scon(j) a(l,imode)=kpar(l,j)*b(l,imode) enddo enddo else if (iabs == 1) then do imode=1,ncon(j)*nspecies !cfpp$ select (vector) do l=1,ldb*scon(j) a(l,imode)=abs(kpar(l,j))*b(l,imode) enddo enddo else if (iabs == 2) then do imode=1,ncon(j)*nspecies !cfpp$ select (vector) do l=1,ldb*scon(j) a(l,imode)=-ii*abs(kpar(l,j))*b(l,imode) enddo enddo else if (iabs == 13) then do imode=1,ncon(j)*nspecies !cfpp$ select (vector) do l=1,ldb*scon(j) a(l,imode)=b(l,imode)/(1+(sgpar*kpar(l,j)*s_par)**4) enddo enddo endif isign=1 ! ! if number of transforms too small, use cfft ! if (ncon(j)*nspecies < 2) then do l=1,ldb*scon(j) c(l)=a(l,1) enddo call cfft(isign,ldb*scon(j),scale,c,1,d,1,ttablekp,ntablekp,work,nwork) do l=1,ldb*scon(j) b(l,1)=d(l) enddo else call xmcfft(isign,ldb*scon(j),ncon(j)*nspecies,scale,a,1,& size(a, 1),b,1,size(b, 1),ttablekp,ntablekp,work,nwork) endif if(iabs < 10) then do k=1,ncon(j)*scon(j) do i=1,nspecies m=mmcon(j,k) n=nncon(j,k) imode=mcon(m,n)+(i-1)*ncon(j) do l=1,ldb f(l,m,n,i)=b(l+ldb*(conpos(m,n)-1),imode)*sgpar enddo enddo enddo else do k=1,ncon(j)*scon(j) do i=1,nspecies m=mmcon(j,k) n=nncon(j,k) imode=mcon(m,n)+(i-1)*ncon(j) do l=1,ldb f(l,m,n,i)=b(l+ldb*(conpos(m,n)-1),imode) enddo enddo enddo endif if (ncycle == 1) then do k=1,ntablekp tablekp(j,k)=ttablekp(k) enddo endif enddo deallocate (ttablekp) else ! ! non-connected way for iperiod=0,1,2 ! scale=1.0/sqrt(float(ldb)) nwork=4*ld*md*nd*nspecies ntablekp=100+2*size(f, 1) do i=1,nspecies do m=1,md do n=1,nd imode=n+(m-1)*nd+(i-1)*md*nd do l=1,ldb a(l,imode)=f(l,m,n,i) enddo enddo enddo enddo isign=-1 call xmcfft(isign,ldb,md*nd*nspecies,scale,a,1,size(a, 1),b,1,size(b,1), & tablekp,ntablekp,work,nwork) if (iabs == -1) then do imode=1,md*nd*nspecies do l=1,ldb a(l,imode)=ii*kpar(l,1)*b(l,imode) enddo enddo else if (iabs == 0) then do imode=1,md*nd*nspecies do l=1,ldb a(l,imode)=kpar(l,1)*b(l,imode) enddo enddo else if (iabs == 1) then do imode=1,md*nd*nspecies do l=1,ldb a(l,imode)=abs(kpar(l,1))*b(l,imode) enddo enddo else if (iabs == 2) then do imode=1,md*nd*nspecies do l=1,ldb a(l,imode)=-ii*abs(kpar(l,1))*b(l,imode) enddo enddo else if (iabs == 13) then do imode=1,md*nd*nspecies do l=1,ldb a(l,imode)=b(l,imode)/(1+(sgpar*kpar(l,1)*s_par)**4) enddo enddo endif isign=1 call xmcfft(isign,ldb,md*nd*nspecies,scale,a,1,size(a, 1),b,1,size(b, 1), & tablekp,ntablekp,work,nwork) if(iabs < 10) then do i=1,nspecies do m=1,md do n=1,nd imode=n+(m-1)*nd+(i-1)*md*nd do l=1,ldb f(l,m,n,i)=b(l,imode)*sgpar enddo enddo enddo enddo else do i=1,nspecies do m=1,md do n=1,nd imode=n+(m-1)*nd+(i-1)*md*nd do l=1,ldb f(l,m,n,i)=b(l,imode) enddo enddo enddo enddo endif endif deallocate (work) end subroutine kparfft