module itg_data ! ! (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. ! implicit none ! ! major physics and control parameters ! integer nspecies ! number of species integer igeo ! switch for general geometry integer iflow ! switch for treatment of m=0 modes integer icrit ! switch for L_Tcrit, D_ML searches integer idt ! switch for adjustable time step integer jrad ! flux surface label (for multiple calls to itgc) real dtadj ! dt -> 1.05*dtadj*dt or dt -> dt/dtadj real dt0 ! time step (in units of Ln/cs) integer lin ! 1= linear run, 0=nonlinear run. integer :: iperiod ! 1 -> periodic in theta b.c.'s, 0 -> 0 b.c.'s ! 2 -> extensions, 3 -> connected FFTs integer :: iseed, n0 integer n_p ! number of particles per cloud integer n_zpic ! number of theta grid points in PIC code real bmax ! bmax of PIC code real chi_int ! int(chi dt) logical binary ! true -> binary i/o for big files logical ascii ! true -> ascii i/o logical cdf ! true -> netcdf i/o real stochf ! stochasticity parameter for poisson logical :: linlay = .false. ! true if layout = 'linear' character layout*6 ! determines layout of k_r character runname*80 ! Name of this run. Associated i/o files ! will be named RUNNAME.in, etc. integer lrunname ! length of runname character note*80 ! note about the run. integer lhistory ! phi @ r(lhistory) saved for plots integer :: mfr, nfr, nfreq, ntotal, ntimf real :: etae ! eta_e real :: rmime ! m_i / m_e real, dimension(:), allocatable :: time, dt ! ! All in the circular concentric small aspect ratio limit unless igeo=1. ! ! omega_gb: omega_gradB vs. theta (toroidal terms) ! omega_kap: omega_kappa vs. theta (toroidal terms) ! omegade: omega_de vs. theta (electron terms) ! b_mag: |B| vs. theta ! bmaginv: 1/B vs. theta ! bgrad: 1/B dB/dl vs. theta real, dimension(:,:,:,:), allocatable :: omega_gb, omega_kap real, dimension(:,:,:), allocatable :: omegade real, dimension(:), allocatable :: b_mag, bmaginv, bgrad ! b_mag is the magnitude of the magnetic field B normalized to B_0 ! gpar is the coefficient of d/dz ( = Ln/qR if igeo=0) ! ! Finite beta ! real beta_e ! 4 Pi ne0 Te0/B0^2 real nuei ! nu_ei * L_n/v_ti * m_e/m_i integer inuei ! flag to select collison model ! 0: (default) simplest Lorentz, -nuei (u_e-u_i) in ! electron momentum eqn only real meovmi ! m_e / m_i used for electron Landau damping in ! EM case. default is zero. integer semi_imp ! ! semi_imp is a parameter for the type of semi-implicit method used ! in the electron moment equations ! 0=no semi-implicit, no predictor corrector in apar eqn ! 1=no semi-implicit, predictor corrector in apar eqn ! 2=weak semi-implicit, predictor corrector in apar eqn ! 3=normal semi-implicit, predictor corrector in apar eqn ! 4=strong semi-implicit, predictor corrector in apar eqn real a0 ! goes with semi_imp real, dimension(:,:,:), allocatable :: omega_gb_e, omega_kap_e real :: cflx, cfly, vmax, rdiff_0 real :: vy, vz, shr, qsf, epsn, eps, epse, alpha, gamma_0, eps_over_q real, dimension(:,:), allocatable :: rdiff real :: etai, nuii, nueeff, rmu1, tol real :: tiovte, cs, etaipar, s_par, s_perp, pfac integer :: igradon, iphi00, iflr, nstp, ifilter, ihdf, nprnt, nread, & ninterv, navgflag, ncycle, ne, nfreqcnt, nperpmom, nparmom, & nemom, iodd, ntrace, iexp, inlpm, stepnlps integer ikx ! Obsolete. Used to control pseudo-spectral vs. ! finite-difference d/dx in nonlinear terms. real :: pmag, dens, avgflag, vxmax, vymax ! ! Geometrical coefficients ! real :: gpar real, dimension(:), allocatable :: gbd, gbd0, cvd, cvd0, jacobian, gd2, gd21, gd22, gradrho real, dimension(10) :: rho, tau, charge, n_I, rmass, eta, eta_par, vt, Ln, nu_ii real, dimension(:,:,:), allocatable :: etak, etak_par ! Closure coefficients: real :: rkkin, rkkperp, rkkperp2, & nuar, nuai, nubr, nubi, nucr, nuci, nudr, nudi, & nu1r, nu1i, nu2r, nu2i, nu3r, nu3i, nu4r, nu4i, & nu5r, nu5i, nu6r, nu6i, nu7r, nu7i, nu8r, nu8i, & nu9r, nu9i, nu10r, nu10i, nu11r, nu11i, nu12r, nu12i, & nu13r, nu13i, nu14r, nu14i, zDq1, zBq1, zEq1 ! When doing debug runs, it is sometimes useful to force the plots ! to have the same date/time label as a previous run, so that one ! can verify results by typing "diff run1.m run2.m". In order to ! do this, set debug_plotlabel to the previous date/time label: character debug_plotlabel*30 contains subroutine input use mp, only: broadcast, proc0 use command_line use gryffin_grid, only: ld, nffty, nfftx, md, nd, x0, xp, y0, z0 use constants, only: pi ! ! Read RUNNAME.in input file, and set some default input values. ! integer i, ierr integer icount character command*80 integer :: movieon ! included for backwards compatibility. Does nothing. integer :: mhi ! included for backwards compatibility. Does nothing. integer :: mlow ! included for backwards compatibility. Does nothing. namelist/wdat/ld,nffty,nfftx,md,nd,dt0,nstp,nprnt,nfreq,nread, & lin,x0,xp,y0,z0,pmag,shr,qsf,epsn,eps,alpha,etai,dens,tiovte, & nueeff,nuii,rmu1,etaipar,vy,vz,dtadj,note,epse,iflow,pfac, & igradon,iphi00,iflr,ihdf,nspecies,etae,rmime, & ninterv,mfr,nfr,nperpmom,nparmom,nemom,iperiod,n0, & iodd,s_par,s_perp,ntrace,ifilter,iexp,ikx, & idt,inlpm,n_p,n_zpic,bmax, & nuar,nuai,nubr,nubi,nucr,nuci,nudr,nudi, & nu1r,nu1i,nu2r,nu2i,nu3r,nu3i,nu4r,nu4i,nu5r,nu5i,nu6r,nu6i, & nu7r,nu7i,nu8r,nu8i,nu9r,nu9i,nu10r,nu10i,nu11r,nu11i, & nu12r,nu12i,nu13r,nu13i,nu14r,nu14i,iseed,stochf, & tau,rmass,charge,eta,eta_par,Ln,n_I,jrad,binary, ascii, cdf, & icrit,gamma_0,rdiff_0,tol,igeo, debug_plotlabel, layout, movieon, & beta_e, semi_imp, nuei, inuei, meovmi, mlow, mhi layout='linear' note='ITG is a nonlinear gyrofluid turbulence simulation code.' gamma_0=0.1 ! set smaller for more accurate L_Tcrit estimate binary=.false. ascii=.false. cdf=.true. tiovte = 1.0 pfac = 0. rdiff_0=0.05 beta_e = 0.0 semi_imp = 0 nuei = 0. inuei = 0 meovmi = 0. icrit=0 igeo=0 jrad=1 n_p=500 n_zpic=32 bmax=1. tol=0.01 nffty=0 nfftx=0 md=0 nd=0 iflow=1 ninterv=10 nfreq=-2 xp=1.e10 x0=1.e10 stochf=0.0 dens=1.0 idt=1 dtadj=0.1 iperiod=0 ! default is old b.c.'s n0=10 epse=0.0 ! r/R for electrons, 0=adiabatic electrons only. igradon=1 s_par=0. ! particle shaping factor s_perp=0. ! particle shaping factor iexp=2 ! used in moment filter ikx=2 inlpm=1 ifilter=0 iphi00=2 iflr=8 rmu1=-1000.0 ! initialization for all diff values iodd=0 ! even and odd modes initialized by default ntrace=0 nperpmom=1 nparmom=3 nemom=3 etaipar=-100. ! arrays: eta = -100. eta_par = -100. tau = -100. rmass = -100. Ln = -100. charge = -100. Ln(1)=1.0 n_I(1)=1.0 lin=1 nspecies=1 rmime=1. etae=0. alpha=0. vy=0.0 vz=0.0 ! set closure coefficients ! ! rkkin is the coefficient for the Hammett-Perkins ! Landau damping model. ! ! rkkin = chi_1 * sqrt(2) * v_ti / c_s ! = 2 * (2/pi*Ti/Te)**0.5, ! ! rkkin=2.*(2./pi)**0.5 rkkperp=(2./pi)**0.5 rkkperp2=(pi/(2.))**0.5 ! ! zDq1 and zBq1 are the coefficients for the 4-pole Landau damping model ! from Greg's PRL: ! ! zDq1 = D1 * sqrt(2) * Sqrt(TiovTe) ! D1 = 2 sqrt(pi) / (3 pi -8) ! zBq1 = (3 + 2*B1)*TiovTe ! zEq1 = (3 + 2*B1) ! B1 = (32 - 9 pi) / (6 pi - 16) zDq1 = 2.*sqrt(pi)/(3.0*pi-8.0)*sqrt(2.) zBq1 = (3.+2.*(32.-9.*pi)/(6.*pi-16.)) zEq1 = 3.+2.*(32.-9.*pi)/(6.*pi-16.) ! Toroidal terms nu1r=-1.e10 ! values which will trigger default nu's. nucr=-1.e10 nudr=-1.e10 ! ! nu11-14 no longer used ! nu11r= 0. nu11i= 0. nu12r= 0. nu12i= 0. nu13r= 0. nu13i= 0. nu14r= 0. nu14i= 0. ! NLPM defaults ! lambda1=1.0 ! lambda2=1.6 ! nu1=1.6 ! nu2=1.3 ihdf=1 ! 1->hdf movie, 0->arrays for particle code. iseed=1 debug_plotlabel=' ' ! determine the name of the input file: if(proc0) then icount=iargc() if(icount >= 1) then call cl_getarg(1,runname,lrunname,ierr) else write(6,*) 'What is the name of the *.in input file?' write(6,*) '(without the .in suffix):' read(5,*) runname lrunname=index(runname,' ')-1 endif write(6,*) 'Running itg for ', runname(1:lrunname) ! Run a shell script to strip out tab characters and comments, which ! some compilers can't handle, from the namelist input file: command='./gfbin/nmlstrip '//runname(1:lrunname)//'.in '//char(0) call ishell(command) ! ! read namelist ! open(unit=12,file=runname(1:lrunname)//'.ins',status='old') read(12,wdat) close(unit=12) if(layout == 'linear') linlay = .true. if(ascii .and. binary) then write(*,*) 'Forcing ascii output' binary = .false. endif if(ascii .and. cdf) then write(*,*) 'Forcing netcdf output' ascii = .false. endif if(cdf .and. binary) then write(*,*) 'Forcing netcdf output' binary = .false. endif if(.not. cdf .and. .not. ascii .and. .not. binary) then write(*,*) 'Choosing NetCDF I/O' cdf = .true. endif endif call broadcast(lrunname) call broadcast(runname) ! Check to see if nspecies > 10. if(nspecies > 10) then if(proc0) write(*,*) 'nspecies > 10 requires recoding the defaults in ' if(proc0) write(*,*) 'the file in.f90. Stopping.' stop endif call broadcast(ld); call broadcast(nffty); call broadcast(nfftx); call broadcast(md) call broadcast(nd); call broadcast(dt0); call broadcast(nstp); call broadcast(nprnt) call broadcast(nfreq); call broadcast(nread); call broadcast(lin); call broadcast(x0) call broadcast(xp); call broadcast(y0); call broadcast(z0); call broadcast(pmag) call broadcast(shr); call broadcast(qsf); call broadcast(epsn); call broadcast(eps) call broadcast(alpha); call broadcast(etai); call broadcast(dens); call broadcast(tiovte) call broadcast(nueeff); call broadcast(nuii); call broadcast(rmu1); call broadcast(etaipar) call broadcast(vy); call broadcast(vz); call broadcast(dtadj); call broadcast(pfac) call broadcast(epse); call broadcast(iflow); call broadcast(igradon); call broadcast(iphi00) call broadcast(iflr); call broadcast(ihdf) call broadcast(nspecies); call broadcast(etae); call broadcast(rmime); call broadcast(ninterv) call broadcast(mfr); call broadcast(nfr); call broadcast(nperpmom); call broadcast(nparmom) call broadcast(nemom); call broadcast(iperiod); call broadcast(n0); call broadcast(iodd) call broadcast(s_par); call broadcast(s_perp); call broadcast(ntrace); call broadcast(ifilter) call broadcast(iexp); call broadcast(ikx); call broadcast(idt); call broadcast(inlpm) call broadcast(n_p); call broadcast(n_zpic); call broadcast(bmax); call broadcast(nuar) call broadcast(nuai); call broadcast(nubr); call broadcast(nubi); call broadcast(nucr) call broadcast(nuci); call broadcast(nudr); call broadcast(nudi); call broadcast(nu1r) call broadcast(nu1i); call broadcast(nu2r); call broadcast(nu2i); call broadcast(nu3r) call broadcast(nu3i); call broadcast(nu4r); call broadcast(nu4i); call broadcast(nu5r) call broadcast(nu5i); call broadcast(nu6r); call broadcast(nu6i); call broadcast(nu7r) call broadcast(nu7i); call broadcast(nu8r); call broadcast(nu8i); call broadcast(nu9r) call broadcast(nu9i); call broadcast(nu10r); call broadcast(nu10i); call broadcast(nu11r) call broadcast(nu11i); call broadcast(nu12r); call broadcast(nu12i); call broadcast(nu13r) call broadcast(nu13i); call broadcast(nu14r); call broadcast(nu14i); call broadcast(iseed) call broadcast(stochf); call broadcast(tau); call broadcast(rmass); call broadcast(charge) call broadcast(eta); call broadcast(eta_par); call broadcast(Ln); call broadcast(n_I) call broadcast(jrad); call broadcast(binary); call broadcast(icrit); call broadcast(gamma_0) call broadcast(rdiff_0); call broadcast(tol); call broadcast(igeo); call broadcast(linlay) call broadcast(ascii); call broadcast(cdf); call broadcast(beta_e); call broadcast(semi_imp) call broadcast(nuei); call broadcast(inuei); call broadcast(meovmi) if(proc0) open(unit=9,file=runname(1:lrunname)//'.out') ! Open output file ! ! Choose multiple of pi for x0 if x0<0: ! if(x0 == 1.e10) then x0=pi*2.5 else if(x0 < 0) x0=pi*abs(x0) endif if(xp == 1.e10) then xp=pi*2.5 else if(xp < 0) xp=pi*abs(xp) endif ! ! If searching for critical gradient, force linear run only: ! if(icrit /= 0 .and. lin == 0) then if(proc0) write(*,*) 'Forcing lin=1 because icrit /= 0' lin=1 endif if(icrit == 2 .and. rmu1 == -1000.) rmu1=0.01 ! Set default value of etaipar if(etaipar == -100.) etaipar=etai ! ! Set self-consistent alpha for finite beta, if alpha = -1. ! correct only for single ion species PS 12/97 ! if (alpha == -1.) & alpha=2.*qsf*qsf*(beta_e/epsn)*(1.+etae+tiovte*(1.+etai)) ! ! By definition these quantities are fixed: ! tau(1)=1.0 rmass(1)=1.0 ! Setup the different rho's, etc. for multiple species: do i=1,nspecies if(eta(i) == -100) eta(i)=etai if(eta_par(i) == -100.) eta_par(i)=eta(i) if(charge(i) == -100.) charge(i)=1. if(Ln(i) == -100.) Ln(i)=1. if(tau(i) == -100.) tau(i)=1. if(rmass(i) == -100.) rmass(i)=1. ! ! Derived quantities: ! rho(i)=sqrt(tau(i)*rmass(i))/charge(i) vt(i)=sqrt(tau(i)/rmass(i)) enddo cs = sqrt(1./tiovte) nu_ii(1)=nuii if(nspecies >= 2) then do i=2,nspecies nu_ii(i)=nuii*charge(i)**4*n_I(i)*vt(i)/(n_I(1)*tau(i)**2) enddo endif ! ! check for errors or problems in the namelist input: if(nstp/nprnt < 2) then if(proc0) write(*,*) 'Warning: Reducing nprnt !' nprnt=nstp/2 endif if(nfreq == -2) nfreq = nprnt ntimf = nfreq nfreqcnt = 0 ! if(iperiod >= 0 .and. iperiod <= 2) then ! indomz=ndomz ! inconz=nconz ! if(indomz /= 1 .or. inconz /= 1) then ! write(*,*) 'Error: ndomz and nconz parameters in itg_data.f90' ! write(*,*) ' MUST be set to 1 for iperiod = 0, 1, or 2' ! call aborter(6,' aborter called from input.f') ! endif ! endif ! ! If closure coefficients aren't set in the namelist, then ! use the values from M.A.Beer Ph.D. thesis. ! if(nu1r == -1.e10) then if(proc0) write(*,*) 'setting default closure coefficients from Beer''s thesis' if (nperpmom == 2.and.nparmom == 4) then nu1r=2.019 nu1i=-1.620 nu2r=.433 nu2i= 1.018 nu3r=-.256 nu3i=1.487 nu4r=-.070 nu4i=-1.382 nu5r=-8.927 nu5i=12.649 nu6r= 8.094 nu6i= 12.638 nu7r= 13.720 nu7i= 5.139 nu8r= 3.368 nu8i= -8.110 nu9r= 1.974 nu9i= -1.984 nu10r= 8.269 nu10i= 2.060 elseif (nperpmom == 1.and.nparmom == 3) then nu1r=1.232 nu1i=.437 nu2r=-.912 nu2i=.362 nu3r=-1.164 nu3i=.294 nu4r=.478 nu4i=-1.926 nu5r=.515 nu5i=-.958 else call aborter(6,'only set up for 3+1 or 4+2 moments') endif endif if(nemom == 3) then if(nucr == -1.e10) then nuar=.290 nuai=-.071 nubr=-1.102 nubi=-.689 nucr=.817 nuci=1.774 endif elseif(nemom == 4) then if(nudr == -1.e10) then nuar=-.038 nuai=.073 nubr=.657 nubi=-.060 nucr=-1.522 nuci=-1.085 nudr=.905 nudi=2.073 endif else call aborter(6,'only set up for 3 or 4 electron moments') endif stepnlps=0 if (proc0 .and. beta_e > 0) then write(*,*) 'beta_e=',beta_e write(*,*) 'semi_imp=',semi_imp write(*,*) 'cs=',cs write(*,*) 'alpha=',alpha write(*,*) 'nuei=',nuei write(*,*) 'inuei=',inuei endif end subroutine input subroutine alloc_itg use gryffin_layouts use gryffin_grid, only: ld, nd, kd, kdpass ! ! Allocate arrays for gryffin. ! allocate ( omega_gb(ld, m_low:m_alloc, n_low:n_alloc, nspecies), & omega_kap(ld, m_low:m_alloc, n_low:n_alloc, nspecies)) omega_gb = 0.; omega_kap = 0. if (beta_e > 0.) then allocate (omega_gb_e(ld, m_low:m_alloc, n_low:n_alloc)) allocate (omega_kap_e(ld, m_low:m_alloc, n_low:n_alloc)) omega_gb_e = 0.; omega_kap_e = 0. endif allocate (time(m_low:m_alloc), dt(m_low:m_alloc)) time = 0.; dt = 0. allocate (b_mag(ld), bmaginv(ld), bgrad(ld), gbd(ld), gbd0(ld), cvd(ld), & cvd0(ld), jacobian(ld), gd2(ld), gd21(ld), gd22(ld), gradrho(ld)) b_mag = 0.; bmaginv = 0.; bgrad = 0.; gbd = 0.; gbd0 = 0.; cvd = 0. cvd= 0.; jacobian = 0.; gd2 = 0.; gd21 = 0.; gd22 = 0.; gradrho = 0. end subroutine alloc_itg subroutine init_itg use gryffin_layouts use gryffin_grid, only: mrr, ld, md, nd, y0, z0 use constants, only: pi integer i if(lin == 1) igradon=0 avgflag=0.0 navgflag=0 ncycle = 0 ne=0 lhistory=float(ld)*.6+1 ntotal = nstp/nfreq ! ! Set up diffusivities: ! if(rmu1 == -1.) ifilter=0 if(icrit == 2) then ifilter=1 rmu1=max(rmu1,0.01) endif allocate (rdiff(m_low:m_alloc, n_low:n_alloc)) rdiff = 0. if(ifilter > 0) then rdiff=rmu1 else rmu1=float(n_p) !so it will appear in old output endif ! ! Find cfl multipliers: ! ! cflx=float(ldb)/x0 ! ! This line gives trouble if nd=1, but we don't run nonlinearly with nd=1, ! so just replace it with something similar: ! ! cflx=shr*z0*float(nd)/(float(nd/2)*2*pi*y0) if(nd == 1) then cflx=shr*z0/(2.*pi*y0) else cflx=shr*z0*float(nd)/(float(nd/2)*2*pi*y0) endif cfly=float(2*(md-1))/(2.*pi*y0) ! ! Allow for k_perp dependent etas for linear estimates: ! allocate (etak(md, nd, nspecies), etak_par(md, nd, nspecies)) do i=1,nspecies etak(:,:,i)=eta(i) etak_par(:,:,i)=eta_par(i) enddo end subroutine init_itg subroutine geo ! ! (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. ! ! Initialize qinp, shat, b_mag, rkperp2, omegad, bgrad for general geometry ! from a run of eiktest. use mp use gryffin_layouts use fitpack_routines use gryffin_grid, only: ld, r, mrr, md, y0, r0n, r0, meq, & rkx, rky, rkperp2, rkx2, rky2 integer :: ntgrid real, dimension(:), allocatable :: theta, bmag, gbdrift, gbdrift0, & gds2, gds21, gds22, cvdrift, cvdrift0, gradpar, logb, tgradrho, dum real :: dpsidrho, drhodpsin, rmaj, bi character*40 char integer :: i, ntg, l, m, n ! ! read in data ! if(proc0) then open(unit=21,file=runname(1:lrunname)//'.geo',status='old') read(21,*) char read(21,*) ntgrid, i, i, drhodpsin, rmaj, shr, bi endif call broadcast(ntgrid) call broadcast(drhodpsin) call broadcast(rmaj) call broadcast(shr) call broadcast(bi) dpsidrho = 1./drhodpsin ntg = 2*ntgrid+1 allocate(theta(ntg), bmag(ntg), gbdrift(ntg), gbdrift0(ntg), & gds2(ntg), gds21(ntg), gds22(ntg), cvdrift(ntg), & cvdrift0(ntg), gradpar(ntg), logb(ntg), tgradrho(ntg)) if(proc0) then read(21,*) char do i=1,ntg read(21,*) gbdrift(i), gradpar(i), tgradrho(i), theta(i) gbdrift(i) = gbdrift(i)*rmaj/4. gradpar(i) = gradpar(i)*rmaj*epsn tgradrho(i) = tgradrho(i)/rmaj/epsn enddo endif call broadcast(gbdrift) call broadcast(gradpar) call broadcast(tgradrho) call broadcast(theta) if(proc0) then read(21,*) char do i=1,ntg read(21,*) cvdrift(i), gds2(i), bmag(i) cvdrift(i) = cvdrift(i)*rmaj/4. logb(i) = log(bmag(i)) enddo endif call broadcast(cvdrift) call broadcast(gds2) call broadcast(bmag) call broadcast(logb) if(proc0) then read(21,*) char do i=1,ntg read(21,*) gds21(i), gds22(i) enddo endif call broadcast(gds21) call broadcast(gds22) if(proc0) then read(21,*) char do i=1,ntg read(21,*) cvdrift0(i), gbdrift0(i) cvdrift0(i) = cvdrift0(i)*rmaj/4. gbdrift0(i) = gbdrift0(i)*rmaj/4. enddo endif call broadcast(cvdrift0) call broadcast(gbdrift0) if(proc0) close(21) allocate (dum(ld)) call inter_d_cspl(ntg,theta,logb,ld,r,dum,bgrad) deallocate (dum) call inter_cspl(ntg,theta,bmag,ld,r,b_mag) call inter_cspl(ntg,theta,gbdrift, ld,r,gbd ) call inter_cspl(ntg,theta,gbdrift0,ld,r,gbd0) call inter_cspl(ntg,theta,cvdrift, ld,r,cvd ) call inter_cspl(ntg,theta,cvdrift0,ld,r,cvd0) call inter_cspl(ntg,theta,gds2, ld,r,gd2 ) call inter_cspl(ntg,theta,gds21,ld,r,gd21) call inter_cspl(ntg,theta,gds22,ld,r,gd22) call inter_cspl(ntg,theta,tgradrho,ld,r,gradrho) gpar = gradpar(1) do l=1,ld bgrad(l)=bgrad(l)*gpar*epsn enddo ! make omegad and rkperp2: !cfpp$ skip R do i=1,nspecies do n = n_low, n_high do m = m_low, m_high do l=1,ld rkx(l,m,n)=(mrr(m,n)/y0)*sqrt(gd22(l))/b_mag(l) ! ! The 1/B above needs to be handled consistently throughout ! the code. ! omega_gb(l,m,n,i)=epsn*rky(m,n)*rho(i)*vt(i) & *(gbd(l)+r0(m,n)*gbd0(l)) omega_kap(l,m,n,i)=epsn*rky(m,n)*rho(i)*vt(i) & *(cvd(l)+r0(m,n)*cvd0(l)) rkperp2(l,m,n)=rky2(m,n) & *(gd2(l)+2.*r0(m,n)*gd21(l)+r0(m,n)**2*gd22(l)) ! ! special cases for k_theta=0, k_r=const. ! if ((m == meq).and.(md /= 1) .and. idx_local(m_lo, meq, n)) then rkx(l,m,n) = -1./y0*r0n(n)*shr/b_mag(l) omega_gb(l,m,n,i)=epsn*rho(i)*vt(i)/y0*r0n(n)*gbd0(l) omega_kap(l,m,n,i)=epsn*rho(i)*vt(i)/y0*r0n(n)*cvd0(l) rkperp2(l,m,n)=(r0n(n)/y0)**2*gd22(l) endif if (md == 1 .and. mrr(m,n) == 0) then rkx(l,m,n)=1./y0 omega_gb(l,m,n,i)=epsn*rho(i)*vt(i)/y0*gbd0(l) omega_kap(l,m,n,i)=epsn*rho(i)*vt(i)/y0*cvd0(l) rkperp2(l,m,n)=1./y0**2 endif ! ! should be krho ! rkx2(l,m,n)=rkx(l,m,n)**2 rkperp2(l,m,n)=rkperp2(l,m,n)/b_mag(l)**2 enddo enddo enddo enddo !cfpp$ noskip R ! ! Define the Jacobian of the coordinate transformation: ! (The abs() is probably overkill.) ! jacobian = abs(dpsidrho/(gpar*b_mag)) bmaginv = 1./b_mag eps_over_q = dpsidrho/bi end subroutine geo subroutine s_alpha use gryffin_layouts use gryffin_grid, only: r, mrr, y0, r0, rkx, rkx2, rkperp2, md, ld, & meq, neq, r0n, rky2, rky integer :: i, l, m, n gpar=epsn/qsf eps_over_q = eps/qsf b_mag = 1./(1.+eps*cos(r)) bmaginv = 1./b_mag bgrad = gpar*eps*sin(r)*b_mag ! ! Backwards compatibility: ! do l=1,ld gradrho(l)=1.0 enddo ! ! define kx grid: ! ! The current version (6.0.1.1) of Cray's fpp has a bug in the following ! nested loops, unless vectorization is turned off: !fpp$ skip R do n = n_low, n_high do m = m_low, m_high ! ! Need switch for m=0 modes. ! rkx(:,m,n)=(mrr(m,n)/y0)*((r(:)-r0(m,n))*shr-alpha*(sin(r(:))-sin(r0(m,n)))) if (m == meq .and. md /= 1 .and. idx_local(m_lo, meq, n)) & rkx(:,m,n)=1./y0*(-r0n(n))*shr if (md == 1 .and. mrr(m,n) == 0) rkx(:,m,n)=1./y0 ! ! local if qsf < 0 ! if (qsf < 0) rkx(:,m,n)=0.0 ! ! rkx(:,m,n)=rkx(:,m,n)*bmaginv(:) ! This will be done properly sometime... ! rkx2(:,m,n)=rkx(:,m,n)**2 ! rkperp2(:,m,n)=rky2(m,n)+rkx2(:,m,n) ! ! should be (k_perp rho)**2 ! rkperp2(:,m,n)=(rky2(m,n)+rkx2(:,m,n))*bmaginv(:)**2 enddo enddo ! !fpp$ noskip R ! ! omega_gb(l) is the grad B drift profile (for toroidal terms) ! omega_kap(l) is the curvature drift profile (for toroidal terms) ! !fpp$ skip R do i=1,nspecies do n = n_low, n_high do m = m_low, m_high omega_gb(:,m,n,i)=0.5*epsn*rky(m,n)*rho(i)*vt(i)* & ( cos(r(:))+(shr*(r(:)-r0(m,n)) & -alpha*(sin(r(:))-sin(r0(m,n))))*sin(r(:))) omega_kap(:,m,n,i) = omega_gb(:,m,n,i) ! ! special case for k_theta=0, k_r=const. ! if ((m == meq).or.((md == 1).and.mrr(m,n) == 0)) then omega_gb(:,m,n,i)=0.5*epsn*rho(i)*vt(i)*rkx(:,m,n)*sin(r(:)) omega_kap(:,m,n,i) = omega_gb(:,m,n,i) if(iflow == 0) then omega_gb(:,m,n,i)=0. omega_kap(:,m,n,i)=0. endif endif ! ! local if qsf<0 ! if (qsf < 0) then omega_gb(:,m,n,i)=epsn*rky(m,n)*rho(i)*vt(i) omega_kap(:,m,n,i) = omega_gb(:,m,n,i) endif enddo enddo enddo !fpp$ noskip R ! ! electron terms in finite beta case ! !fpp$ skip R if (beta_e > 0) then do n = n_low, n_high do m = m_low, m_high omega_gb_e(:,m,n)=0.5*epsn*rky(m,n)*(-1.0/tiovte)* & ( cos(r(:))+(shr*(r(:)-r0(m,n)) & -alpha*(sin(r(:))-sin(r0(m,n))))*sin(r(:))) omega_kap_e(:,m,n) = omega_gb_e(:,m,n) ! ! special case for k_theta=0, k_r=const. ! if ((m == meq).or.((md == 1).and.mrr(m,n) == 0)) then omega_gb_e(:,m,n)=0.5*epsn*(-1.0/tiovte)*rkx(:,m,n)*sin(r(:)) omega_kap_e(:,m,n) = omega_gb_e(:,m,n) if(iflow == 0) then omega_gb_e(:,m,n)=0. omega_kap_e(:,m,n)=0. endif endif ! ! local if qsf<0 ! if (qsf < 0) then omega_gb_e(:,m,n)=epsn*rky(m,n)*(-1.0/tiovte) omega_kap_e(:,m,n) = omega_gb_e(:,m,n) endif enddo enddo endif !fpp$ noskip R ! ! Define the Jacobian of the coordinate transformation: ! (The abs() is probably overkill.) ! jacobian = abs(1./(gpar*b_mag)) ! Debug: ! ! i=1 ! do n=1,nd ! do m=1,md ! do l=1,ld ! write(*,2000) l,m,n,i,omegad(l,m,n,i) ! write(*,2001) l,bgrad(l),epsn*eps/qsf*sin(r(l))*b_mag(l),b_mag(l) ! write(*,2002) l,m,n,i,rkperp2(l,m,n) ! write(*,2003) l,gpar,gradrho(l),r(l) ! write(*,*) l,gbd0(l)+cvd0(l),-shr*sin(r(l)) ! write(*,*) l,gbd(l)+cvd(l),cos(r(l))+shr*r(l)*sin(r(l)) ! write(*,*) ! enddo ! write(*,*) ! enddo ! enddo 2000 format('wd(',i2,',',i2,',',i2,',',i2,')= ',f13.6) 2001 format('grad ln(B)(',i2,')= ',f13.6,' ',f13.6,' ',f13.6) 2002 format('kperp**2(',i2,',',i2,',',i2,',',i2,')= ',f13.6) 2003 format('gradpar(',i2,')= ',f13.6,' |grad.rho|= ',f13.6,' theta= ',f13.6) end subroutine s_alpha end module itg_data