subroutine flrinit(init_crit) c c initializes some arrays used in flr, poisson, timestep, and nlpsc c and kparfft c c (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, c and Gregory W. Hammett. ALL RIGHTS RESERVED. c implicit none include 'itg.par' include 'itg.cmn' real rkp2, rgammas, tmu, wgt, zfilter real s18ceE, s18cfE, znorm integer i, j, l, m, n, ifail, lmax, init_crit c c beginning of flr and poisson initializations c flr c******************************************************************** c if(init_crit.eq.0) then do i=1,nspecies do n=1,nd do m=1,md do l=1,ld-1 gnwgt(l,m,n,i)= . (g0wgt(l,m,n,i)+etak(m,n,i)*flrwgt(l,m,n,i)) gtpwgt(l,m,n,i)=etak_par(m,n,i)*g0wgt(l,m,n,i) gtwgt(l,m,n,i)= . etak(m,n,i)*(g0wgt(l,m,n,i)+flrwgt2(l,m,n,i)) . +flrwgt(l,m,n,i) grwgt(l,m,n,i)= . flrwgt(l,m,n,i)+etak(m,n,i)*(g0wgt(l,m,n,i) . +flrwgt2(l,m,n,i)+flrwgt(l,m,n,i)) enddo enddo enddo enddo return endif c******************************************************************** c initialize some arrays only once per run: c ifail=0 do n=1,nd do m=1,md do l=1,ld-1 pfilter2(l,m,n)=0. enddo enddo enddo c c BD: could reverse i,n loops for better multitasking performance c do i=1,nspecies do n=1,nd do m=1,md do l=1,ld-1 rkp2=rkperp2(l,m,n)*rho(i)**2 rgammas=s18cfE(rkp2,ifail)/s18ceE(rkp2,ifail) c c Use "new" Pade approximation c if (iflr.eq.1) then g0wgt(l,m,n,i)=1./(1.+rkp2/2.) flrwgt(l,m,n,i)=-rkp2/(2.+rkp2) flrwgt2(l,m,n,i)=-rkp2/(1.+rkp2/2.)**2 c c Use "new" FLR model c else if(iflr.eq.2.or.iflr.eq.3.or.iflr.eq.8.or. * iflr.eq.9) then g0wgt(l,m,n,i)=sqrt(s18ceE(rkp2,ifail)) flrwgt(l,m,n,i)=-rkp2/2.*(1.-rgammas) flrwgt2(l,m,n,i)=-rkp2/2.*(1.+(1.-rgammas) * *(1.-rkp2/2.*(3.+rgammas))) c c Use exp(-b/2) model c else if(iflr.eq.4) then g0wgt(l,m,n,i)=exp(-rkp2/2.) flrwgt(l,m,n,i)=-rkp2/2. flrwgt2(l,m,n,i)=-rkp2*(1.-rkp2/4.) c c Try Taylor series c else if(iflr.eq.5) then flrwgt(l,m,n,i)=-rkp2/2. flrwgt2(l,m,n,i)=-rkp2/2. g0wgt(l,m,n,i)=1.0 c c Like Ron, but Linsker treatment different c else if(iflr.eq.6) then g0wgt(l,m,n,i)=s18ceE(rkp2,ifail) flrwgt(l,m,n,i)=-rkp2*(1.-rgammas) flrwgt2(l,m,n,i)=-rkp2*(2.-rgammas) * +2.*rkp2**2*(1-rgammas) c c Full FLR, but no Linsker effect c else if(iflr.eq.7) then g0wgt(l,m,n,i)=s18ceE(rkp2,ifail) flrwgt(l,m,n,i)=-rkp2*(1.-rgammas) flrwgt2(l,m,n,i)=-rkp2*(1.-rgammas) endif flrwgt3(l,m,n,i)=flrwgt(l,m,n,i)**2 gnwgt(l,m,n,i)=(1.+etak(m,n,i)*flrwgt(l,m,n,i))*g0wgt(l,m,n,i) gtpwgt(l,m,n,i)=etak_par(m,n,i)*g0wgt(l,m,n,i) gtwgt(l,m,n,i)=etak(m,n,i)*(1.+flrwgt2(l,m,n,i))+flrwgt(l,m,n,i) gtwgt(l,m,n,i)=gtwgt(l,m,n,i)*g0wgt(l,m,n,i) grwgt(l,m,n,i)=(flrwgt(l,m,n,i)+etak(m,n,i)*(1.+flrwgt2(l,m,n,i) * +flrwgt(l,m,n,i)))*g0wgt(l,m,n,i) c c toroidal FLR terms, keep old forms for backward compatability c c Pade in Poisson with new toroidal FLR terms, c suggested forms in M.A.Beer Ph.D. thesis c if (iflr.eq.8) then gnwgtd(l,m,n,i)=(2.+flrwgt(l,m,n,i))*g0wgt(l,m,n,i)/tau(i) gtpwgtd(l,m,n,i)=2.*g0wgt(l,m,n,i)/tau(i) gtwgtd(l,m,n,i)=(1.+2*flrwgt(l,m,n,i) & +flrwgt2(l,m,n,i))*g0wgt(l,m,n,i)/tau(i) gqwgtb(l,m,n,i)=(flrwgt2(l,m,n,i)-flrwgt(l,m,n,i)) & *g0wgt(l,m,n,i) c c Non-O(b) accurate closure in thesis c else if(iflr.eq.9) then gnwgtd(l,m,n,i)=(2.+3.*flrwgt(l,m,n,i) & -flrwgt2(l,m,n,i))*g0wgt(l,m,n,i)/tau(i) gtpwgtd(l,m,n,i)=2.*g0wgt(l,m,n,i)/tau(i) gtwgtd(l,m,n,i)=(1.+flrwgt(l,m,n,i) & +4.*flrwgt2(l,m,n,i))*g0wgt(l,m,n,i)/tau(i) gqwgtb(l,m,n,i)=0.0 ! only tested this model with gpar B=0 c c Old forms of toroidal terms, for backward comp. only c else gnwgtd(l,m,n,i)=2.*g0wgt(l,m,n,i)/tau(i) gtpwgtd(l,m,n,i)=2.*g0wgt(l,m,n,i)/tau(i) gtwgtd(l,m,n,i)=(1.+flrwgt(l,m,n,i)*3.) & *g0wgt(l,m,n,i)/tau(i) gqwgtb(l,m,n,i)=0.0 endif gnwgtd(l,m,n,i)=gnwgtd(l,m,n,i)*charge(i) gtpwgtd(l,m,n,i)=gtpwgtd(l,m,n,i)*charge(i) gtwgtd(l,m,n,i)=gtwgtd(l,m,n,i)*charge(i) flrwgt(l,m,n,i)=flrwgt(l,m,n,i)*g0wgt(l,m,n,i) flrwgt2(l,m,n,i)=flrwgt2(l,m,n,i)*g0wgt(l,m,n,i) flrwgt3(l,m,n,i)=flrwgt3(l,m,n,i)*g0wgt(l,m,n,i) enddo enddo enddo enddo c c end of initialization section c******************************************************************** c c poisson c******************************************************************** c initialize some arrays only once per run: c c precalculating the g0wgt and g2wgt arrays in FLR, and the c pfilter array in POISSON, reduces the CPU time spent in these c subroutines by about 25%. c do i=1,nspecies do n=1,nd do m=1,md do l=1,ld-1 rkp2=rkperp2(l,m,n)*rho(i)**2 c c Use "new" Pade approximation if iflr=1: c if (iflr.eq.1) then wgt=(1.-1./(1.+rkp2)) g0wgtp(l,m,n,i)=2./(2.+rkp2) gewgt(l,m,n,i)=2./(2.+rkp2) twgt(l,m,n,i)=-rkp2/(2.+rkp2) nwgt(l,m,n,i)=1. c c or use "new" FLR model for iflr=2 or 3: c else if(iflr.eq.2.or.iflr.eq.3) then wgt=(1.-s18ceE(rkp2,ifail)) rgammas=s18cfE(rkp2,ifail)/s18ceE(rkp2,ifail) tmu=1.-rkp2*(1.-rkp2/2.)+rkp2/2.* * rgammas*(1.-rkp2*rgammas) twgt(l,m,n,i)=-rkp2/2.*(1.-rgammas) nwgt(l,m,n,i)=1.-rkp2*(1.-rkp2/4.)+rkp2* * rgammas/2.*(1.+rkp2*(1.-3.*rgammas/2.)) gewgt(l,m,n,i)=sqrt(s18ceE(rkp2,ifail)) g0wgtp(l,m,n,i)=gewgt(l,m,n,i)/tmu c c = exp(-b/2) for iflr=4: c else if(iflr.eq.4) then wgt=(1.-s18ceE(rkp2,ifail)) g0wgtp(l,m,n,i)=exp(-rkp2/2.) gewgt(l,m,n,i)=exp(-rkp2/2.) twgt(l,m,n,i)=-rkp2/2. nwgt(l,m,n,i)=1. c c or use Taylor series: c else if(iflr.eq.5) then wgt=rkp2 g0wgtp(l,m,n,i)=1. gewgt(l,m,n,i)=1. twgt(l,m,n,i)=-rkp2/2. nwgt(l,m,n,i)=1. c c Like Ron, but Linsker treatment different c else if(iflr.eq.6) then wgt=(1.-s18ceE(rkp2,ifail)) twgt(l,m,n,i)=0.0 nwgt(l,m,n,i)=1. g0wgtp(l,m,n,i)=1. c c Full FLR, but no Linsker effect c else if(iflr.eq.7) then wgt=(1.-s18ceE(rkp2,ifail)) twgt(l,m,n,i)=0.0 nwgt(l,m,n,i)=1. g0wgtp(l,m,n,i)=1. c c Pade in Poisson with new toroidal FLR terms c else if(iflr.eq.8.or.iflr.eq.9) then wgt=(1.-s18ceE(rkp2,ifail)) g0wgtp(l,m,n,i)=1. gewgt(l,m,n,i)=1. twgt(l,m,n,i)=-rkp2/2./(1.+rkp2/2.)**2 nwgt(l,m,n,i)=1./(1.+rkp2/2.) endif c g0wgtp(l,m,n,i)=g0wgtp(l,m,n,i)*charge(i)*n_I(i) pfilter2(l,m,n)=pfilter2(l,m,n) . +wgt*n_I(i)*charge(i)**2/tau(i) pfilter3(l,m,n,i)=wgt*n_I(i)*charge(i)**2/tau(i) c BD: Shouldn't these rky2 and rkx2 factors have an additional factor of c 1/B**2??? c c GWH: c Handle finite-size particle shapes/filtering. [Should check exactly c at what stage the particle simulations put in their phi filters, i.e., c before or after calculating ? Also, do they filter in real space c or in field-line coordinates? (Doesn't matter in the shearless c limit). Do different species have different sizes?] c c Use Dimit's style of filter (the parallel filter is elsewhere): zfilter=exp(-s_perp**4*rky2(m,n)**2)/ . (1+s_perp**4*rkx2(l,m,n)**2) c c BD: Question: shouldn't this rkperp2 definition have a factor of 1/B**2? c ?? c Standard filter style: c rkperp2=rky2(m,n)+rkx2(l,m,n) c zfilter=exp(-s_perp**2*rkperp2) twgt(l,m,n,i)=twgt(l,m,n,i)*zfilter nwgt(l,m,n,i)=nwgt(l,m,n,i)*zfilter c Put particle shapes/filtering in the coefficients used to find n_i_bar: enddo enddo enddo enddo c do n=1,nd do m=1,md do l=1,ld-1 pfilter(l,m,n)=1.0/(tiovte+pfilter2(l,m,n)) enddo enddo enddo cfpp$ skip R if(iphi00.eq.2 .and. (lin.eq.1 .or. shr.eq.0.)) then do n=1,nd do m=1,md denom0(m,n)=0. if ((rky2(m,n).eq.0.0).and.(rkx2(1,m,n).ne.0.0)) then do l=l_left,l_right-1 denom0(m,n)=denom0(m,n) . +pfilter(l,m,n)*pfilter2(l,m,n)*jacobian(l) c added dl/B factor for flux surface average enddo endif enddo enddo else do n=2,nd denom0(1,n)=0. do l=l_left,l_right-1 denom0(1,n)=denom0(1,n) . +pfilter(l,1,n)*pfilter2(l,1,n)*jacobian(l) enddo enddo endif cfpp$ noskip R c end of initialization section c********************************************************************* c c for kparfft c do j=1,ndom lmax=ldb*scon(j) do l=1,lmax if (l.le.lmax/2) then kpar(l,j)=(l-1)*2*pi/(x0*scon(j)*2.) else kpar(l,j)=(l-1-lmax)*2*pi/(x0*scon(j)*2.) endif c c local if qsf<0 c if (qsf.lt.0 .and. igeo.eq.0) kpar(l,j)=shr enddo enddo return end