subroutine flrinit c c initializes some arrays used in flr, poisson, timestep, and nlpsc c and kparfft c implicit none include 'itg.par' include 'itg.cmn' real rkperp2, rgammas,tmu,wgt real s18ceE,s18cfE integer l,m,n,ifail,index,lmax c c c define ky grid c ifail=0 do n=1,nd do m=1,md rky(m,n)=(mrr(m,n)/y0) ! ky rky2(m,n)=(mrr(m,n)/y0)**2 ! ky**2 enddo enddo c no de-aliasing c c define kx grid: c c do n=1,nd do m=1,md c c Need switch for m=0 modes. c do l=1,ld-1 rkx(l,m,n)=(mrr(m,n)/y0)*((r(l)-r0(m,n))*shr & -alpha*(sin(r(l))-sin(r0(m,n)))) if (m.eq.meq) then rkx(l,m,n)=mrr(m+1,n)/y0*(-r0(m+1,n))*shr endif if (md.eq.1) then rkx(l,m,n)=1./y0 endif c c local if qsf < 0 c if (qsf.lt.0) rkx(l,m,n)=0.0 rkx(l,m,n)=rkx(l,m,n)*bmaginv(l) rkx2(l,m,n)=rkx(l,m,n)**2 enddo enddo enddo c c beginning of flr and poisson initializations c c flr c******************************************************************** c initialize some arrays only once per run: c ifail=0 c index=0 do 18 l=1,ld-1 do 18 n=1,nd do 18 m=1,md index=index+1 rkperp2=rky2(m,n)+rkx2(l,m,n) kxwgt(index)=rkx(l,m,n) c Use "new" Pade approximation if (iflr.eq.1) then g0wgt(index)=1./(1.+rkperp2/2.) flrwgt(index)=-rkperp2/(2.+rkperp2) flrwgt3(index)=flrwgt(index)**2 flrwgt2(index)=-rkperp2/(1.+rkperp2/2.)**2 gnwgt(index)=(1.+etai*flrwgt(index))*g0wgt(index) gtpwgt(index)=etaipar*g0wgt(index) gtwgt(index)=(etai*(1.+flrwgt2(index))+flrwgt(index)) * *g0wgt(index) grwgt(index)=(flrwgt(index)+etai*(1.+flrwgt2(index) * +flrwgt(index)))*g0wgt(index) flrwgt(index)=flrwgt(index)*g0wgt(index) flrwgt2(index)=flrwgt2(index)*g0wgt(index) flrwgt3(index)=flrwgt3(index)*g0wgt(index) c c Use "new" FLR model c else if(iflr.eq.2.or.iflr.eq.3) then rgammas=s18cfE(rkperp2,ifail)/s18ceE(rkperp2,ifail) g0wgt(index)=sqrt(s18ceE(rkperp2,ifail)) flrwgt(index)=-rkperp2/2.*(1.-rgammas) flrwgt3(index)=flrwgt(index)**2 flrwgt2(index)=-rkperp2/2.*(1.+(1.-rgammas)*(1.-rkperp2/2.* * (3.+rgammas))) gnwgt(index)=(1.+etai*flrwgt(index))*g0wgt(index) gtpwgt(index)=etaipar*g0wgt(index) gtwgt(index)=etai*(1.+flrwgt2(index))+flrwgt(index) gtwgt(index)=gtwgt(index)*g0wgt(index) grwgt(index)=(flrwgt(index)+etai*(1.+flrwgt2(index) * +flrwgt(index)))*g0wgt(index) flrwgt(index)=flrwgt(index)*g0wgt(index) flrwgt2(index)=flrwgt2(index)*g0wgt(index) flrwgt3(index)=flrwgt3(index)*g0wgt(index) c c Use exp(-b/2) model c else if(iflr.eq.4) then g0wgt(index)=exp(-rkperp2/2.) flrwgt(index)=-rkperp2/2. flrwgt3(index)=flrwgt(index)**2 flrwgt2(index)=-rkperp2*(1.-rkperp2/4.) gnwgt(index)=(1.+etai*flrwgt(index))*g0wgt(index) gtwgt(index)=(etai*(1.+flrwgt2(index))+flrwgt(index)) * *g0wgt(index) gtpwgt(index)=etaipar*g0wgt(index) grwgt(index)=(flrwgt(index)+etai*(1.+flrwgt2(index) * +flrwgt(index)))*g0wgt(index) flrwgt(index)=flrwgt(index)*g0wgt(index) flrwgt2(index)=flrwgt2(index)*g0wgt(index) flrwgt3(index)=flrwgt3(index)*g0wgt(index) c c Try Taylor series c else if(iflr.eq.5) then flrwgt(index)=-rkperp2/2. flrwgt3(index)=flrwgt(index)**2 flrwgt2(index)=-rkperp2/2. g0wgt(index)=1.0 gnwgt(index)=(1.-rkperp2/2.-etai*rkperp2/2.) gtpwgt(index)=etaipar gtwgt(index)=-rkperp2/2.+etai*(1-3.*rkperp2/2.) c c Like Ron, but Linsker treatment different c else if(iflr.eq.6) then rgammas=s18cfE(rkperp2,ifail)/s18ceE(rkperp2,ifail) g0wgt(index)=s18ceE(rkperp2,ifail) flrwgt(index)=-rkperp2*(1.-rgammas) flrwgt3(index)=flrwgt(index)**2 flrwgt2(index)=-rkperp2*(2.-rgammas)+2.*rkperp2**2*(1-rgammas) gnwgt(index)=(1.+etai*flrwgt(index))*g0wgt(index) gtpwgt(index)=etaipar*g0wgt(index) gtwgt(index)=etai*(1.+flrwgt2(index))+flrwgt(index) gtwgt(index)=gtwgt(index)*g0wgt(index) c grwgt(index)=(flrwgt(index)+etai*(1.+flrwgt2(index) c * +flrwgt(index)))*g0wgt(index) grwgt(index)=(0.5-1.75*rkperp2+1.375*rkperp2*rgammas & +0.75*rkperp2**2*(1-rgammas) )*g0wgt(index) flrwgt(index)=(1.-rkperp2/2.*(1-rgammas))*g0wgt(index) flrwgt2(index)=flrwgt2(index)*g0wgt(index) flrwgt3(index)=flrwgt3(index)*g0wgt(index) c c Full FLR, but no Linsker effect c else if(iflr.eq.7) then rgammas=s18cfE(rkperp2,ifail)/s18ceE(rkperp2,ifail) g0wgt(index)=s18ceE(rkperp2,ifail) flrwgt(index)=-rkperp2*(1.-rgammas) flrwgt3(index)=flrwgt(index)**2 flrwgt2(index)=-rkperp2*(1.-rgammas) gnwgt(index)=(1.+etai*flrwgt(index))*g0wgt(index) gtpwgt(index)=etaipar*g0wgt(index) gtwgt(index)=etai*(1.+flrwgt2(index))+flrwgt(index) gtwgt(index)=gtwgt(index)*g0wgt(index) grwgt(index)=(flrwgt(index)+etai*(1.+flrwgt2(index) * +flrwgt(index)))*g0wgt(index) flrwgt(index)=flrwgt(index)*g0wgt(index) flrwgt2(index)=flrwgt2(index)*g0wgt(index) flrwgt3(index)=flrwgt3(index)*g0wgt(index) endif 18 continue 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 ifail=0 index=0 do 19 l=1,ld-1 do 19 n=1,nd do 19 m=1,md index=index+1 rkperp2=rky2(m,n)+rkx2(l,m,n) c c Use "new" Pade approximation if iflr=1: c if (iflr.eq.1) then if (nspecies.eq.2) then wgt=(1.-1./(1.+rkperp2)) else wgt=tiovte+(1.-1./(1.+rkperp2)) endif g0wgtp(index)=2./(2.+rkperp2) gewgt(index)=2./(2.+rkperp2) twgt(index)=-rkperp2/(2.+rkperp2) nwgt(index)=1. c c or use "new" FLR model for iflr=2 or 3: c else if(iflr.eq.2.or.iflr.eq.3) then if(nspecies.eq.2) then wgt=(1.-s18ceE(rkperp2,ifail)) else wgt=tiovte+(1.-s18ceE(rkperp2,ifail)) endif rgammas=s18cfE(rkperp2,ifail)/s18ceE(rkperp2,ifail) c tmu=1.-rkperp2*(1.-rkperp2/2.)+rkperp2/2.* * rgammas*(1.-rkperp2*rgammas) twgt(index)=-rkperp2/2.*(1.-rgammas) nwgt(index)=1.-rkperp2*(1.-rkperp2/4.)+rkperp2* * rgammas/2.*(1.+rkperp2*(1.-3.*rgammas/2.)) gewgt(index)=sqrt(s18ceE(rkperp2,ifail)) g0wgtp(index)=gewgt(index)/tmu c c = exp(-b/2) for iflr=4: c else if(iflr.eq.4) then if(nspecies.eq.2) then wgt=(1.-s18ceE(rkperp2,ifail)) else wgt=tiovte+(1.-s18ceE(rkperp2,ifail)) endif g0wgtp(index)=exp(-rkperp2/2.) gewgt(index)=exp(-rkperp2/2.) twgt(index)=-rkperp2/2. nwgt(index)=1. c c or use Taylor series: c else if(iflr.eq.5) then if(nspecies.eq.2) then wgt=rkperp2 else wgt=tiovte+rkperp2 endif g0wgtp(index)=1. gewgt(index)=1. twgt(index)=-rkperp2/2. nwgt(index)=1. c c Like Ron, but Linsker treatment different c else if(iflr.eq.6) then if(nspecies.eq.2) then wgt=(1.-s18ceE(rkperp2,ifail)) else wgt=tiovte+(1.-s18ceE(rkperp2,ifail)) endif c twgt(index)=0.0 nwgt(index)=1. g0wgtp(index)=1. c c Full FLR, but no Linsker effect c else if(iflr.eq.7) then if(nspecies.eq.2) then wgt=(1.-s18ceE(rkperp2,ifail)) else wgt=tiovte+(1.-s18ceE(rkperp2,ifail)) endif c twgt(index)=0.0 nwgt(index)=1. g0wgtp(index)=1. endif c pfilter2(index)=wgt c c If k_perp^2 = 0 and nonadiabatic electrons are being used, c then make Phi for that k_perp = 0. c if(nspecies.eq.2) then if(wgt.gt.1.e-33) then pfilter(index)=1./wgt else pfilter(index)=0. endif else pfilter(index)=1./wgt endif c c Filter out high kx's: c if(beta.ne.0.) * pfilter(index)=pfilter(index)*exp(-beta**2*rkperp2) 19 continue c end of initialization section c********************************************************************* c c for kparfft c lmax=ld-1 do l=1,lmax if (l.le.lmax/2) then kpar(l)=(l-1)*2*pi/(x0*2.*qsf/epsn) else kpar(l)=(l-1-lmax)*2*pi/(x0*2.*qsf/epsn) endif c c local if qsf<0 c if (qsf.lt.0) kpar(l)=shr enddo return end