module nl ! ! (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, ! 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 itg_data implicit none private public :: nlpsc, init_nl, courant integer :: idum real, dimension(:,:), allocatable, save :: kxa2, kya2, kxb2, kyb2 real, dimension(:,:), allocatable :: km real, dimension(:), allocatable :: work real, dimension(:), allocatable :: worky real, dimension(:), allocatable :: tablex ! mcfft trig array real, dimension(:), allocatable :: tabley ! mcfft trig array integer :: ntablex contains subroutine nlpsc(isave, a, b, nl_term, lmax) 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(ld,md,nd), kya(ld,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(ld,md,nd), kyb(ld,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. ! 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, nwork, ntablex real :: scale 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 interface subroutine csfftm(isign, n, lot, scale, x, ldx, & y, ldy, table, work, isys) integer, intent(in) :: isign, n, lot, ldx, ldy integer, intent(in) :: isys real, intent(in) :: scale complex, dimension(0:ldx-1,0:lot-1), intent(in) :: x real, dimension(0:ldy-1,0:lot-1), intent(out) :: y real, dimension(100+2*n) :: table real, dimension((2*n+4)*lot) :: work end subroutine csfftm end interface nwork=4*lmax*md*nalias ntablex=100+2*nalias scale=1. isign=1 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 allocate(work(nwork)) call xmcfft(isign,nalias,md*lmax,scale,ka,1,nalias,ka,1, & nalias,tablex,ntablex,work,nwork) deallocate(work) ! ! 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 allocate (worky((2*malias+4)*lmax*nalias)) 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) 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,nwork,ntablex integer l,m,n,index,imode 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 interface subroutine scfftm(isign, n, lot, scale, x, ldx, & y, ldy, table, work, isys) integer, intent(in) :: isign, n, lot, ldx, ldy integer, intent(in) :: isys real, intent(in) :: scale complex, dimension(0:ldy-1,0:lot-1), intent(in) :: y real, dimension(0:ldx-1,0:lot-1), intent(out) :: x real, dimension(100+2*n) :: table real, dimension((2*n+4)*lot) :: work end subroutine scfftm end interface scale = 1./float(malias) isign = -1 allocate (worky((2*malias+4)*lmax*nalias)) 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 nwork=4*lmax*md*nalias scale=1./float(nalias) ntablex=100+2*nalias isign=-1 allocate(work(nwork)) call xmcfft(isign,nalias,md*lmax,scale,aky,1,nalias,aky,1,nalias, & tablex,ntablex,work,nwork) deallocate(work) 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 end subroutine k_space subroutine init_nl ! initialize trig tables for y FFT's: idum=0 allocate (tabley(100+2*malias)) call csfftm(0,malias,1,1.0,idum,idum,idum,idum,tabley,idum,0) ! initialize trig tables for x FFT's: ntablex=100+2*nalias allocate (tablex(ntablex)) call xmcfft(0,nalias,1,1.0,0.0,idum,idum,0.0,idum,idum,tablex,ntablex,1.0,idum) end subroutine init_nl subroutine mod_alloc(lmax) 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 courant use flr_terms, only: phi_u real :: vmax1 integer :: i, m ! ! ******************************************************************* ! Check Courant condition: ! vmax1=0. do i=1,nspecies call maxv(phi_u(:,:,:,i),vmax1) ! ! With our FLR models, the nonlinear FLR terms are always smaller than ! the main term, so the next two lines are unnecessary. This may not be ! true for such models as the Taylor series FLR with nonlinear FLR. ! ! call maxv(phi_flr1(:,:,:,i),vmax1) ! call maxv(phi_flr2(:,:,:,i),vmax1) enddo vmax1=max(dtadj/dt0,vmax1) do m=1,md dt(m)=min(dt0,dtadj/vmax1) enddo end subroutine courant subroutine maxv(b,max_vel) complex, dimension(:,:,:) :: b complex, dimension(:,:,:), pointer :: kxb, kyb real max_vel integer l,m,n,imode,mr allocate (kxb(ld,md,nd), kyb(ld,md,nd)) 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) deallocate (kxb) call real_space(kyb,kyb2,ldb) deallocate (kyb) 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 end module nl