module gryffin_grid ! ! (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. ! ! ! This module contains routines and variables related to the computational grid. ! implicit none integer ld ! # of parallel meshpoints integer ldb ! (# of parallel mespoints) - 1 integer kd ! # of kappa grid points (pitch angle) trapped integer kdpass ! # of kappa grid points (pitch angle) +passing integer nffty ! # of y grid points, with dealiasing integer nfftx ! # of x grid points, with dealiasing integer nmin ! fourier modes range from nmin to nmax. integer nmax ! nmin=nmax=0 for single helicity integer nd ! total # of n's, nd=nmax-nmin+1 integer :: md, nddp, nddm, malias, nalias, madd, nadd, l_left, l_right, meq, neq real :: x0, xp, y0, z0, rootn, rootm real, dimension(:), allocatable :: r, dr! grid for the radial (x) coordinate real, dimension(:), allocatable :: kap ! grid for kappa (pitch angle) coordinate real, dimension(:,:), allocatable :: r0 real, dimension(:), allocatable :: r0n ! Many versions ago, a wedge of modes was kept; mmin, mmax specified the wedge. ! nrr and mrr contain the actual n and m numbers for a given index value. ! integer, dimension(:), allocatable :: mmin, mmax integer, dimension(:), allocatable :: nrr integer, dimension(:,:), allocatable :: mrr, mrr_all ! ! FOURIER TRANSFORM CONVENTIONS FOR THE ITG CODE ! ! ?? The description below is for the old slab version of the code, ! and needs updating for the toroidal version: ! ! The code solves fluid equations using a spectral representation for ! the poloidal (y) and toroidal (z) directions, and a second-order ! finite difference method for the radial (x) direction. All ! quantities can be expanded in Fourier harmonics as: ! ! wfield(x,y,z) = Sum_mn w(x,m,n) exp( + i m y / y0 - i n z / z0) ! ! where w(x,m,n) are complex coefficients. ! ! The (-m,-n) mode is linearly dependent on the (m,n) mode, so ! we can eliminate half of the modes. We choose to keep just ! the modes with m >= 0, although n <0, =0, and >0 are all be kept. ! Caution: though the m>0 modes are all independent, the m=0 modes ! are not. ! ! All primary variables are for guiding centers, normalized to the ! equilibrium value for each species, i.e. n(j)=(L_ne/rho_i) n_j1/n_j0, ! except the potential, which is Phi=(L_ne/rho_i) e Phi_1/T_i0, where ! T_i0 is the main species equilibrium temperature, rho_i=v_ti/Omega_i, ! and v_ti=sqrt(T_i0/m_i) (E.g., 'density' is the guiding center density, ! which is different from the particle density due to FLR effects.) ! ! Time is normalized to L_ne/v_ti. ! ! density = guiding center density n ! u_par = parallel velocity ! t_par = parallel temperature ! q_par = parallel flux of parallel heat ! t_perp = perpendicular temperature ! q_perp = parallel flux of perpendicular heat ! potential = phi (electrostatic potential) ! ! The "home" basis for most fields is (theta,ky,kx). For example, ! density(lz,mz,nz,nspecz) is the guiding center density vs. ! l=1,ldb along-the-field-line grid index ! m=1,md ky (poloidal mode) index ! n=1,nd kx (radial mode) index ! nsp=1,nspecies species index real, dimension(:,:), allocatable :: rky, rky2 ! ky, ky**2 for each mode real, dimension(:,:,:), allocatable :: rkx, rkx2, rkperp2 ! kr, kr**2, kperp**2 for each mode real, dimension(:,:,:), allocatable :: rkperp20,rkperp21,rkperp22 contains subroutine alloc_grid use gryffin_layouts allocate (r(ld), dr(ld)) r = 0.; dr = 0. allocate (mmin(nd), mmax(nd)) mmin = 0; mmax = 0 allocate (r0(m_low:m_alloc, n_low:n_alloc), r0n(n_low:n_alloc)) r0 = 0.; r0n = 0. allocate (mrr(m_low:m_alloc, n_low:n_alloc), mrr_all(md, nd)) mrr = 0; mrr_all = 0 allocate (rkx(ld, m_low:m_alloc, n_low:n_alloc), & rkx2(ld, m_low:m_alloc, n_low:n_alloc), & rkperp2(ld, m_low:m_alloc, n_low:n_alloc)) allocate (rkperp20(ld, m_low:m_alloc, n_low:n_alloc), & rkperp21(ld, m_low:m_alloc, n_low:n_alloc), & rkperp22(ld, m_low:m_alloc, n_low:n_alloc)) rkx = 0.; rkx2 = 0.; rkperp2 = 0. rkperp20 = 0. ; rkperp21 = 0. ; rkperp22 = 0. ; allocate (rky(m_low:m_alloc, n_low:n_alloc), & rky2(m_low:m_alloc, n_low:n_alloc)) rky = 0.; rky2 = 0. end subroutine alloc_grid subroutine init_grid (epse, iperiod, lin, linlay) use gryffin_layouts use mp, only: proc0, send, receive, barrier real, intent (in) :: epse integer, intent (in) :: iperiod, lin logical, intent (in) :: linlay real dum integer :: l, m, n, mdspec ldb=ld-1 ! ! initialize arrays of m and n modes being used: ! ! nd # of n modes ! md # of m modes for each n ! malias=0 nalias=0 if(nffty /= 0) malias=nffty if(md == 0) md=(nffty-1)/3+1 if(nffty /= 0) nalias=nfftx if(nd == 0) nd=2*((nfftx-1)/3)+1 if(nd == -1) then nmin=0 nmax=1 nddm=1 nddp=1 allocate(nrr(1)) nrr=1 nd=1 goto 413 endif if(mod(nd,2) == 0) nd=nd+1 if(nd == 1) then nmin=0 nmax=0 nddm=1 nddp=1 allocate(nrr(1)) nrr=0 else allocate(nrr(nd)) nddm=nd/2 nddp=nd/2+1 nmin=-nd/2 nmax=nd/2 do n=1,nddp nrr(n)=n-1 enddo do n=1,nddm nrr(n+nddp)=n-nddm-1 enddo endif 413 continue if(epse > 0.0) then ! This criterion for kd needs to be upgraded for the case x0>pi, as does ! much of the trapped electron calculation: kd=ldb/2 ! covers just the trapped particles kdpass=ldb ! covers all particles, including passing particles else ! initialize to a definite value, since these are used to define some ! array sizes irrespectives of the value of epse: kd=0 kdpass=0 endif ! Set up for real-to-complex y FFT's. ! md = # of independent ky modes (with ky>=0): mdspec=2*(md-1) ! = # of ky points (without dealiasing) if(malias == 0) malias=4*(md-1) ! overkill (more efficient to specify nffty) ! if(malias > mffty) call aborter(6, & ! 'error: nffty or malias > mffty (see itg.in and itg_data.f90)') madd=malias-mdspec rootm=sqrt(float(malias)) ! ! Set up for complex-to-complex x FFT's. ! nd = # of independent kx modes (usually odd for symmetry around kx=0) ! nalias = # of x points (with dealiasing): if(nd == 1) nalias=1 if(nalias == 0) nalias=2*(nd-1) ! overkill (more efficient to specify nfftx) ! if(nalias > mfftx) call aborter(6, & ! 'error: nfftx or nalias > mfftx (see itg.in and itg_data.f90)') nadd=nalias-nd rootn=sqrt(float(nalias)) call init_gryffin_layouts (md, nd) call alloc_grid ! r(l) is the grid for the radial coordinate: do l=1,ld r(l)=2.*x0*float(l-1)/ldb-x0 enddo do l=2,ldb dr(l)=r(l+1)-r(l-1) enddo dr(1)=2.*(r(2)-r(1)) dr(ld)=2.*(r(ld)-r(ldb)) ! neq forced to be one and I did away with possibility of lowest mrr being non-zero ! Bill 10.20.98 mmax = md -1 mmin = 0 do n=n_low, n_high do m=m_low, m_high mrr(m,n)=m+mmin(n)-1 enddo enddo ! ! assign mrr_all: ! do m=1, md mrr_all(m,:)=m-1 enddo ! ! define ky grid ! rky = mrr/y0 rky2 = (mrr/y0)**2 ! BD: removed 1.5.99 (temporarily, I hope) ?? ! if(proc0) then ! write(9,wdat) ! write(9,1100) (nrr(n),mmin(n),mmax(n),n=1,nd) ! endif 1100 format(3x,'at n =',i4,5x,'mmin =',i4,3x,'mmax =',i4) ! ! Twisted periodicity boundary conditions parameters: ! if (iperiod == 2) then l_left=float(ld-1)/2.*(1.-xp/x0)+1.+.5 l_right=float(ld-1)/2.*(1.+xp/x0)+1.+.5 if(proc0) write(9,*) 'l_left, l_right: ', l_left, l_right else l_left=1 l_right=ld endif ! ! r0=theta0, mbeer ! if(proc0) write(9,*) 'Theta_0:' do m = m_low,m_high do n = n_low, n_high if (mrr(m,n) == 0.or.(lin == 1 .and. linlay)) then if (nmax /= 0) then r0(m,n)=z0*float(nrr(n))/nmax else r0(m,n)=z0*float(nrr(n)) endif else if (nmax /= 0) then r0(m,n)=z0*float(nrr(n))/nmax/float(mrr(m,n)) else r0(m,n)=z0*float(nrr(n))/float(mrr(m,n)) endif endif ! define theta0 for the k_theta=0 modes: if (m == 1) then if (nmax /= 0) then r0n(n)=z0*float(nrr(n))/nmax else r0n(n)=z0*float(nrr(n)) endif endif enddo enddo ! ! Send all values of r0 to processor zero for output ! do n = m_lo%nl_world, m_lo%nu_world do m = m_lo%ml_world, m_lo%mu_world if(proc0) then if(idx_local(m_lo, m, n)) then dum = r0(m,n) else call receive(dum, proc_id(m_lo, m, n)) endif write(9,*) m,n,' r0= ',dum elseif (idx_local(m_lo, m, n)) then call send(r0(m,n), 0) endif call barrier enddo enddo ! indices of k = 0 mode meq = 1 neq = 1 end subroutine init_grid end module gryffin_grid