module nlps ! ! (c) Copyright 1991 to 1998 by Michael A. Beer, William D. Dorland, ! P. B. Snyder, Q. P. Liu, and Gregory W. Hammett. ALL RIGHTS RESERVED. ! ! "NonLinear PseudoSpectral" BD ! Given the fields a, b, ka, and kb this subroutine calculates the ! Poisson bracket ! {a,b} = da/dx db/dy - da/dy db/dx. ! ! The fields arrive in (x,ky,kz) space and must be transformed ! with the appropriate derivatives into (x,y,z) space. ! use fft_work implicit none private public :: nlpsc, init_nlps, maxv integer :: idum real, dimension(:,:), allocatable, save :: kxa2, kya2, kxb2, kyb2 real, dimension(:,:), allocatable :: km real, dimension(:), allocatable :: workx real, dimension(:), allocatable :: worky real, dimension(:), allocatable :: tablex ! mcfft trig array real, dimension(:), allocatable :: tabley ! mcfft trig array integer :: ntablex, ntabley, nworkx, nworky contains subroutine nlpsc(isave, a, b, nl_term, lmax) use constants, only: ii use gryffin_grid, only: mrr, y0, r0 integer l, m, n, isave, mr, lmax integer, save :: lmax_last integer :: nfirst = 1 complex, dimension(:,:,:) :: a, b, nl_term complex, dimension(:,:,:), allocatable :: kxa, kya, kxb, kyb ! Don't recalculate FFT's unnecessarily! ! ! isave determines which fields are new for this call: ! isave=0 -> both fields are new ! isave=1 -> "a" fields are the same as last time nlps was called ! isave=2 -> "b" fields are the same as last time nlps was called ! if(nfirst == 1) then nfirst = 0 call mod_alloc(lmax) else if(lmax /= lmax_last) then deallocate (kxa2, kya2, kxb2, kyb2, km)!, km1) call mod_alloc(lmax) endif endif lmax_last = lmax ! Calculate the d/dy and d/dx terms: if(isave /= 1) then allocate (kxa(lmax,md,nd), kya(lmax,md,nd)) do n=1,nd do l=1,lmax kya(l,1,n)=0. kxa(l,1,n)=ii/y0*shr*r0(1,n)*a(l,1,n) enddo enddo do n=1,nd do m=2,md mr=mrr(m,n) do l=1,lmax kya(l,m,n)=ii*mr/y0*a(l,m,n) kxa(l,m,n)=ii*mr/y0*shr*r0(m,n)*a(l,m,n) enddo enddo enddo endif if(isave /= 2) then allocate (kxb(lmax,md,nd), kyb(lmax,md,nd)) do n=1,nd do l=1,lmax kyb(l,1,n)=0. kxb(l,1,n)=ii/y0*shr*r0(1,n)*b(l,1,n) enddo enddo do n=1,nd do m=2,md mr=mrr(m,n) do l=1,lmax kyb(l,m,n)=ii*mr/y0*b(l,m,n) kxb(l,m,n)=ii*mr/y0*shr*r0(m,n)*b(l,m,n) enddo enddo enddo endif ! Transform the vectors i*kx*a, i*kx*b, i*ky*b and i*ky*a into x space: if(isave /= 1) then call real_space(kxa, kxa2, lmax) deallocate (kxa) call real_space(kya, kya2, lmax) deallocate (kya) endif if(isave /= 2) then call real_space(kxb,kxb2,lmax) deallocate (kxb) call real_space(kyb,kyb2,lmax) deallocate (kyb) endif ! Perform the multiplications: km = kxa2 * kyb2 - kxb2 * kya2 ! Transform back. call k_space(nl_term, km, lmax) end subroutine nlpsc subroutine real_space(ak, a, lmax) ! ! Transform a variable from k space into real space and dealias. ! use gryffin_grid, only: nddp, nadd, nalias, md, malias complex, dimension(:,:,:) :: ak real, dimension(:,:) :: a integer :: lmax complex, dimension(nalias,md*lmax) :: ka complex, dimension(malias/2+1,nalias*lmax) :: aky_tmp integer :: l, m, n, index, imode, isign real :: scale do n=1,nddp do m=1,md do l=1,lmax imode=l+(m-1)*lmax ka(n,imode)=ak(l,m,n) enddo enddo enddo do n=nadd+nddp+1,nalias do m=1,md do l=1,lmax imode=l+(m-1)*lmax ka(n,imode)=ak(l,m,n-nadd) enddo enddo enddo ! ! de-alias ! do n=nddp+1,nadd+nddp do m=1,md do l=1,lmax imode=l+(m-1)*lmax ka(n,imode)=0.0 enddo enddo enddo scale=1. isign=1 nworkx = ccfftm_work0 + ccfftm_work1*lmax*md*nalias allocate(workx(nworkx)) call mcfft(isign,nalias,md*lmax,scale,ka,1,nalias,ka,1, & nalias,tablex,ntablex,workx,nworkx) deallocate(workx) ! ! Get the ky=0 (m=1) component: ! do n=1,nalias do l=1,lmax imode=l+(n-1)*lmax aky_tmp(1,imode)=ka(n,l) enddo enddo do n=1,nalias do m=2,md do l=1,lmax imode=l+(n-1)*lmax index=l+(m-1)*lmax aky_tmp(m,imode)=ka(n,index)/2. enddo enddo enddo ! ! dealias for high ky: ! aky_tmp(md+1:malias/2+1,:)=0.0 nworky = (csfftm_work1*malias+csfftm_work2)*lmax*nalias allocate (worky(nworky)) call csfftm(isign,malias,lmax*nalias,scale,aky_tmp,malias/2+1, & a,malias,tabley,worky,0) deallocate (worky) end subroutine real_space subroutine k_space(ak, ay, lmax) use gryffin_grid, only: malias, md, nalias, nddp, nadd real, dimension(:,:) :: ay complex, dimension(:,:,:) :: ak integer :: lmax complex, dimension(nalias, md*lmax) :: aky complex, dimension(malias/2+1,nalias*lmax) :: aky_tmp real scale integer isign integer l,m,n,index,imode scale = 1./float(malias) isign = -1 nworky = (csfftm_work1*malias+csfftm_work2)*lmax*nalias allocate (worky(nworky)) call scfftm(isign,malias,lmax*nalias,scale,ay,malias, & aky_tmp,malias/2+1,tabley,worky,0) deallocate (worky) ! ! get ky=0 (m=1) component: ! do n=1,nalias do l=1,lmax imode=l+(n-1)*lmax aky(n,l)=aky_tmp(1,imode) enddo enddo ! ! get ky>0 components (factor of 2 for sin/cos vs. complex conventions). ! do n=1,nalias do m=2,md do l=1,lmax imode=l+(n-1)*lmax index=l+(m-1)*lmax aky(n,index)=aky_tmp(m,imode)*2. enddo enddo enddo isign=-1 scale=1./float(nalias) nworkx = ccfftm_work0 + ccfftm_work1*lmax*md*nalias allocate(workx(nworkx)) call mcfft(isign,nalias,md*lmax,scale,aky,1,nalias,aky,1,nalias, & tablex,ntablex,workx,nworkx) deallocate(workx) do n=1,nddp do m=1,md do l=1,lmax imode=l+(m-1)*lmax ak(l,m,n)=aky(n,imode) enddo enddo enddo do n=nadd+nddp+1,nalias do m=1,md do l=1,lmax imode=l+(m-1)*lmax ak(l,m,n-nadd)=aky(n,imode) enddo enddo enddo ! BD: 1.5.99 ?? ! we need to enforce reality condition here ! since the calls elsewhere have been taken away call reality(ak) end subroutine k_space subroutine init_nlps complex zdum ! initialize trig tables for y FFT's: idum=0 ntabley = csfftm_table0 + csfftm_table1*malias allocate (tabley(ntabley)) call csfftm(0,malias,1,1.0,idum,idum,idum,idum,tabley,idum,0) ! initialize trig tables for x FFT's: ntablex = ccfftm_table0 + ccfftm_table1*nalias allocate (tablex(ntablex)) call mcfft(0,nalias,1,1.0,zdum,idum,idum,zdum,idum,idum,tablex,ntablex,1.0,idum) end subroutine init_nlps subroutine mod_alloc(lmax) use gryffin_grid, only: malias, nalias integer :: lmax allocate (kxa2(malias, nalias*lmax), kya2(malias, nalias*lmax)) allocate (kxb2(malias, nalias*lmax), kyb2(malias, nalias*lmax)) allocate (km(malias, nalias*lmax)) end subroutine mod_alloc subroutine maxv(b,max_vel) use constants, only: ii use gryffin_grid, only: ldb, md, nd, mrr, malias, nalias, y0, r0 use itg_data complex, dimension(:,:,:) :: b complex, dimension(ld,md,nd) :: kxb, kyb real max_vel integer l,m,n,imode,mr do n=1,nd do l=1,ldb kyb(l,1,n)=0. kxb(l,1,n)=ii/y0*shr*r0(1,n)*b(l,1,n) enddo enddo do m=2,md do n=1,nd mr=mrr(m,n) do l=1,ldb kyb(l,m,n)=ii*mr/y0*b(l,m,n) kxb(l,m,n)=ii*mr/y0*shr*r0(m,n)*b(l,m,n) enddo enddo enddo deallocate (kxb2, kyb2) allocate (kxb2(malias, nalias*ldb), kyb2(malias, nalias*ldb)) call real_space(kxb,kxb2,ldb) call real_space(kyb,kyb2,ldb) do m=1,malias do imode=1,nalias*ldb max_vel=max(max(max_vel,abs(kyb2(m,imode))*cflx), & abs(kxb2(m,imode))*cfly) enddo enddo end subroutine maxv subroutine reality (a) complex, dimension(:,:,:) :: a use gryffin_grid, only: nddm, nddp, nadd integer :: m, n m = 1 do n = 1, nddm a(:,m,nddp+nadd+n) = conjg(a(:,m,nddp+1-n) enddo end subroutine reality end module nlps