subroutine prstep(t_start,dt,time) c c Prints out various quantities of interest. c Checks energy conservation, only valid for eps=0. c This calculates the analytic formula c for dW/dt, and compares with W(t+dt)-W(t-dt). c implicit none include 'itg.par' include 'itg.cmn' c real dt(mz),time(mz),phisqold(mz,nz),phirms1 real dz(lz),dz0,znorm, t_start, deltat real dum,dummy(mz,nz),nuur,nuui,den2,v2,tpar2,tp2,vol real dum1k(mz,nz),dum2k(mz,nz) real den2m(mz,nz),v2m(mz,nz),tpar2m(mz,nz),tp2m(mz,nz),pr2m(mz,nz) real d1eta,d2eta,d3eta,d1par,d2par,wupsi real d1r,d3r,d4r,d5r,d6r,d7r,d8r,d9r real d1eta2(mz,nz),d2eta2(mz,nz),d3eta2(mz,nz),d1par2(mz,nz) real d2par2(mz,nz),wupsi2(mz,nz),d7r2(mz,nz),d8r2(mz,nz) real d9r2(mz,nz) real d1r2(mz,nz),d3r2(mz,nz),d4r2(mz,nz),d5r2(mz,nz),d6r2(mz,nz) complex totpr(lz,mz,nz,nspecz),ky_psiv(lz,mz,nz,nspecz) complex ky_phi(lz,mz,nz),ky_phiba(kz,mz,nz) complex nbar(lz,mz,nz,nspecz) complex ky_psiq(lz,mz,nz,nspecz),p_perp(lz,mz,nz,nspecz) complex kyn1(lz,mz,nz,nspecz),kyt1(lz,mz,nz,nspecz) complex kyt2(lz,mz,nz,nspecz),kwn1(lz,mz,nz,nspecz) complex kwt1(lz,mz,nz,nspecz),kwt2(lz,mz,nz,nspecz) complex kpt1(lz,mz,nz,nspecz),kpt2(lz,mz,nz,nspecz) complex ktu1(lz,mz,nz,nspecz),ktt1(lz,mz,nz,nspecz) complex ktt2(lz,mz,nz,nspecz),ktr1(lz,mz,nz,nspecz) complex ktt3(lz,mz,nz,nspecz),ktt4(lz,mz,nz,nspecz) complex bcpot integer ninv(-nz:nz),nshift,n_max,n1,j1 integer i,l,m,n,init,j,j0,mr,nr,mpts,npts real ky,flux1,flux2,x,y,beta1,bz data init /1/ save init c if (nparmom.eq.3) then nuur=nu5r nuui=nu5i else nuur=0.0 nuui=0.0 endif dtold=dt(md) c c new growth rate stuff , mbeer if(init.eq.0) then do n=1,nd do m=1,md phisqold(m,n)=phisq(ne-1,m,n) enddo enddo endif call volsq(potential,potential,phirms1,dummy) do m=1,md do n=1,nd phisq(ne,m,n)=dummy(m,n) enddo enddo phirms(ne)=sqrt(phirms1) do n=1,nd do m=1,md if ((phisqold(m,n).ne.0.) .and. (init.eq.0)) then mgamx(m,n,ne)=log(phisq(ne,m,n)/phisqold(m,n)) & /2./dt(m)/float(nprnt) else mgamx(m,n,ne)=0. endif enddo enddo c c end of growth rate stuff c c Evaluate drive terms again to make dE/dt 2nd order in dt c c c Driving terms c c ********************* c c load arrays to be volume averaged c call flr(potential) do i=1,nspecies do m=1,md do n=1,nd do l=1,ld c slab drive kyn1(l,m,n,i)=ii*rky(m,n)*phi_n(l,m,n,i)/Ln(i) kyt1(l,m,n,i)=ii*rky(m,n)*phi_tperp(l,m,n,i)/Ln(i) kyt2(l,m,n,i)=ii*rky(m,n)*phi_tpar(l,m,n,i)/Ln(i) c toroidal drive kwn1(l,m,n,i)=-ii*omegad(l,m,n,i)*phi_nd(l,m,n,i) kwt1(l,m,n,i)=-2.*ii*omegad(l,m,n,i)* . (nu2i*t_perp(l,m,n,i)+0.5*phi_tpard(l,m,n,i) ) kwt2(l,m,n,i)=-2.*ii*omegad(l,m,n,i)* . ( nu3i*t_par(l,m,n,i)+0.5*phi_tperpd(l,m,n,i)) c parallel damping kpt1(l,m,n,i)=rkkin*t_par(l,m,n,i)*vt(i) kpt2(l,m,n,i)=rkkperp* . (t_perp(l,m,n,i)+phi_qperp(l,m,n,i) . *charge(i)/tau(i))*vt(i) c toroidal damping ktu1(l,m,n,i)=2.*nu5r*abs(omegad(l,m,n,i))*u_par(l,m,n,i) ktt1(l,m,n,i)=2.*nu1r*abs(omegad(l,m,n,i))*t_par(l,m,n,i) ktt2(l,m,n,i)=2.*nu2r*abs(omegad(l,m,n,i))*t_perp(l,m,n,i) ktt3(l,m,n,i)=2.*nu3r*abs(omegad(l,m,n,i))*t_par(l,m,n,i) ktt4(l,m,n,i)=2.*nu4r*abs(omegad(l,m,n,i))*t_perp(l,m,n,i) c transfer ktr1(l,m,n,i)=phi_u(l,m,n,i)*vt(i)*charge(i)/tau(i) enddo enddo enddo enddo call kparfft(kpt1,1) call kparfft(kpt2,1) call kparfft(ktr1,-1) c c Do the volume averages, add to drive c c BD igeo=1 changes, 6/25/95: c I haven't worked through this section carefully. It may be that I c am wrong to replace volsq calls with volflux calls here. Since c this is just a diagnostic, I will come back to it later. c c olddrive=drive c call volflux(density,kyn1,d1eta,d1eta2) c call volflux(t_par,kyt2,d2eta,d2eta2) c call volflux(t_perp,kyt1,d3eta,d3eta2) c c ?? BD 6/25/95: c Not generalized for igeo=1 call volsq(u_par,ktr1,wupsi,wupsi2) c c ?? BD 6/25/95: c Not generalized for igeo=1 call volsq(t_par,kpt1,d1par,d1par2) c c ?? BD 6/25/95: c Not generalized for igeo=1 call volsq(t_perp,kpt2,d2par,d2par2) c call volflux(density,kwn1,d1r,d1r2) c call volflux(t_par,kwt1,d4r,d4r2) c call volflux(t_perp,kwt2,d7r,d7r2) c call volflux(u_par,ktu1,d3r,d3r2) c nu**r call volflux(t_par,ktt1,d5r,d5r2) c nu**r call volflux(t_par,ktt2,d6r,d6r2) c nu**r call volflux(t_perp,ktt4,d9r,d9r2) c nu**r call volflux(t_perp,ktt3,d8r,d8r2) drive=-d1eta-d1r-wupsi-d3r-0.5*(d1par+d2eta+d4r+d5r+d6r) . -d2par-d3eta-d7r-d8r-d9r drt(1,ne)=-d1eta drt(2,ne)=-0.5*d2eta drt(3,ne)=-d3eta drt(4,ne)=-0.5*d1par drt(5,ne)=-d2par drt(6,ne)=-d1r drt(7,ne)=-d3r drt(8,ne)=-d4r/2. drt(9,ne)=-d5r/2. drt(10,ne)=-d6r/2. drt(11,ne)=-d7r drt(12,ne)=-d8r drt(13,ne)=-d9r drt(14,ne)=-wupsi c c ******************* c Check energy conservation for 3+1 c c call volsq(density,density,den2,den2m) call volsq(u_par,u_par,v2,v2m) call volsq(t_par,t_par,tpar2,tpar2m) call volsq(t_perp,t_perp,tp2,tp2m) wenr=0.5*(den2+v2+tpar2/2.+tp2) wenx(ne)=wenr do m=1,md do n=1,nd wenrk(ne,m,n)=0.5*(den2m(m,n)+v2m(m,n)+tpar2m(m,n)/2. . +tp2m(m,n)) dketa(ne,m,n)=-d1eta2(m,n)-0.5*d2eta2(m,n)-d3eta2(m,n) dkpar(ne,m,n)=-0.5*d1par2(m,n)-d2par2(m,n) dktor(ne,m,n)=-d1r2(m,n)-d4r2(m,n)-d7r2(m,n) dktdp(ne,m,n)=-d3r2(m,n)-d5r2(m,n)/2.-d6r2(m,n)/2. . -d8r2(m,n)-d9r2(m,n) wkups(ne,m,n)=-wupsi2(m,n) dktot(ne,m,n)=dketa(ne,m,n)+dkpar(ne,m,n) . +dktor(ne,m,n)+dktdp(ne,m,n)+wkups(ne,m,n) enddo enddo c c Total growth rate c if(wenr.ne.0.) gamx(ne)=(wenr-pwenr)/(dt(md)*wenr) c c **************** c c This is the fastest way to code the new flux. c c********************************************************************* c Calculate the heat flux, using the full pressure c where P = n (T_par/2 + T_perp). Thus, in perturbed quantities, c P_1 = 3/2 n_1 + T_perp1 + 1/2 T_par1 c Define Pressures c if (iflr.ne.6) then do i=1,nspecies do n=1,nd do m=1,md ky=mrr(m,n)/y0 do l=1,ld totpr(l,m,n,i)=(t_perp(l,m,n,i)+ * 0.5*t_par(l,m,n,i)+1.5*density(l,m,n,i)) p_perp(l,m,n,i)=density(l,m,n,i)+t_perp(l,m,n,i) ky_psiv(l,m,n,i)=ii*ky*phi_u(l,m,n,i) ky_psiq(l,m,n,i)=ii*ky*phi_qperp(l,m,n,i) enddo enddo enddo enddo c c Do the volume averages c do i=1,nspecies call volflux(totpr(1,1,1,i),ky_psiv(1,1,1,i),flux1,dummy) do n=1,nd do m=1,md qfluximn(ne,m,n,i)=dummy(m,n)*n_I(i)*tau(i)*tiovte enddo enddo call volflux(p_perp(1,1,1,i),ky_psiq(1,1,1,i),flux2,dummy) do n=1,nd do m=1,md qfluximn(ne,m,n,i)=qfluximn(ne,m,n,i)+dummy(m,n) . *n_I(i)*tau(i)*tiovte enddo enddo wpfx(ne,i)=(flux1+flux2)*n_I(i)*tau(i)*tiovte enddo if(ne.eq.1) then deltat=time(md)-t_start else deltat=time(md)-timo(ne-1) endif chi_int=chi_int-wpfx(ne,1)*deltat/eta(1)+1.5*rdiff(1,1)*deltat c else c c Ron's way c do i=1,nspecies do n=1,nd do m=1,md ky=mrr(m,n)/y0 do l=1,ld totpr(l,m,n,i)=t_perp(l,m,n,i)+ * 0.5*t_par(l,m,n,i)+1.5*density(l,m,n,i) & -1.5*potential(l,m,n)+phi_n(l,m,n,i)/etai & +(1.5-1./etai)*phi_u(l,m,n,i) ky_psiv(l,m,n,i)=ii*ky*potential(l,m,n) enddo enddo enddo enddo c c Do the volume averages c call volflux(totpr,ky_psiv,flux1,dummy) c wpfx(ne,1)=flux1 endif c c Particle and electron heat fluxes c do i=1,nspecies do n=1,nd do m=1,md ky=mrr(m,n)/y0 do l=1,ld nbar(l,m,n,i)=((density(l,m,n,i)*nwgt(l,m,n,i) . +t_perp(l,m,n,i)*twgt(l,m,n,i))*g0wgtp(l,m,n,i) . -pfilter3(l,m,n,i)*potential(l,m,n))/charge(i) ky_psiv(l,m,n,i)=ii*ky*potential(l,m,n) enddo enddo enddo enddo do i=1,nspecies call volflux(nbar(1,1,1,i),ky_psiv(1,1,1,i),flux1,dummy) fluxi(ne,i)=flux1 do m=1,md do n=1,nd fluximn(ne,m,n,i)=dummy(m,n) enddo enddo enddo c do n=1,nd do m=2,md ky=mrr(m,n)/y0 do l=1,kd ky_phiba(l,m,n)=ii*ky*phi_ba(l,m,n) enddo enddo m=1 ky=mrr(m,n)/y0 do l=1,kdpass ky_phiba(l,m,n)=ii*ky*phi_ba(l,m,n) enddo enddo c ?? BD 6/25/95: c Not generalized for igeo=1 (kappa integrals ignored so far in changes) call volkflux(e_density,ky_phiba,flux1,dum1k) call volkflux(phi_ba,ky_phiba,flux2,dum2k) fluxe(ne)=flux1-flux2*tiovte do m=1,md do n=1,nd fluxemn(ne,m,n)=dum1k(m,n)-tiovte*dum2k(m,n) enddo enddo c ?? BD 6/25/95: c Not generalized for igeo=1 call volkflux(e_p,ky_phiba,flux1,dum1k) qfluxe(ne)=1.5*(flux1-flux2*tiovte) do m=1,md do n=1,nd qfluxemn(ne,m,n)=1.5*(dum1k(m,n)-tiovte*dum2k(m,n)) enddo enddo c ******************* c c total energy, Kinetic energy, and pressure vs. mode c (Not really kinetic energy, just v^2 ) c c Define total pressure c call volsq(totpr,totpr,dum,pr2m) do m=1,md do n=1,nd wtif(ne,m,n)=0.5*(den2m(m,n)+v2m(m,n) > +0.5*tpar2m(m,n)+tp2m(m,n)) wkif(ne,m,n)=v2m(m,n) wpif(ne,m,n)=pr2m(m,n) enddo enddo c c Phi* Phi for each m and n at z=0 c l=ldb/2 do m=1,md do n=1,nd psp(ne,m,n)=real(potential(l,m,n))**2+ & aimag(potential(l,m,n))**2 enddo enddo c c Other variables we want to save c timo(ne)=time(md) c c Flux surface averaged poloidal and toroidal flows, for m=0,n=0 mode c if(iperiod.ne.2) then l_left=1 l_right=ld-1 endif dz0=2.*x0/ldb do l=l_left,l_right c backwards compatibility issue: if(igeo.eq.0) then dz(l)=dz0*jacobian(l) else c BD: c factor of jacobian missing for igeo=0, which I am including c to improve the HRW work. c 4.6.98 dz(l)=dz0*jacobian(l) endif enddo utor(ne)=0.0 upol(ne)=0.0 c c ?? BD 6/25/95: c Not generalized for igeo=1 c rkx, bz not upgraded for general geometry yet, so skip this part c if(igeo.eq.0) then bz=1/sqrt(1+eps**2/qsf**2) do i=1,nspecies znorm=0. do l=l_left,l_right utor(ne)=utor(ne)+dz(l)*(real(u_par(l,1,1,i)) & +eps/qsf*bmaginv(l)*bz*rkx(l,1,1)*aimag(phi_u(l,1,1,i))) & *bmaginv(l)**2*bz upol(ne)=upol(ne) & +dz(l)*(real(u_par(l,1,1,i))*eps/qsf & -bmaginv(l)*bz*rkx(l,1,1)*aimag(phi_u(l,1,1,i))) & *bmaginv(l)*bz*bmaginv(l) c BD: This is the original coding: c znorm=znorm+dz(l) c replaced with: znorm=znorm+dz(l)*bmaginv(l) c 4.6.98 BD enddo enddo if(znorm.eq.0.) write(*,*) 'znorm = 0: problem.' utor(ne)=utor(ne)/znorm upol(ne)=upol(ne)/znorm endif ccccccccccc c c flux surface averaged quantities do n=1,nd phi00(ne,n)=0.0 den00(ne,n)=0.0 upar00(ne,n)=0.0 tpar00(ne,n)=0.0 tperp00(ne,n)=0.0 znorm=0. do l=l_left,l_right-1 phi00(ne,n)=phi00(ne,n)+potential(l,1,n)*dz(l) den00(ne,n)=den00(ne,n)+density(l,1,n,1)*dz(l) upar00(ne,n)=upar00(ne,n)+u_par(l,1,n,1)*dz(l) tpar00(ne,n)=tpar00(ne,n)+t_par(l,1,n,1)*dz(l) tperp00(ne,n)=tperp00(ne,n)+t_perp(l,1,n,1)*dz(l) znorm=znorm+dz(l) enddo phi00(ne,n)=phi00(ne,n)/znorm den00(ne,n)=den00(ne,n)/znorm upar00(ne,n)=upar00(ne,n)/znorm tpar00(ne,n)=tpar00(ne,n)/znorm tperp00(ne,n)=tperp00(ne,n)/znorm enddo ccccccccccc c c Correlation functions mpts=4*md npts=4*nd l=ldb/2+1 do j=1,npts cx(ne,j)=0.0 x=2.0*pi*float(j-1-npts/2)/float(npts-1) do m=1,md do n=1,nd nr=nrr(n) mr=mrr(m,n) cx(ne,j)=cx(ne,j)+real(potential(l,m,n)*exp(ii*nr*x)) enddo enddo enddo do j=1,mpts cy(ne,j)=0.0 y=2.0*pi*float(j-1-mpts/2)/float(mpts-1) do m=1,md do n=1,nd nr=nrr(n) mr=mrr(m,n) cy(ne,j)=cy(ne,j)+real(potential(l,m,n)*exp(-ii*mr*y)) enddo enddo enddo c ??BD 9/22/95: c Apparently these integrals need Jacobian factors if igeo=1, c but I am leaving it for later: c j0=ldb/2+1 do j=1,ld cz(ne,j)=0.0 czc(ne,j)=0.0 czcn(ne,j)=0.0 czn(ne,j)=0.0 cznn(ne,j)=0.0 do n=1,nd czn(ne,j)=czn(ne,j) & +real(potential(j,1,n)-phi00(ne,n)) & *real(potential(j0,1,n)-phi00(ne,n)) & +aimag(potential(j,1,n)-phi00(ne,n)) & *aimag(potential(j0,1,n)-phi00(ne,n)) cznn(ne,j)=cznn(ne,j) & +real(potential(j,1,n)-phi00(ne,n)) & *real(potential(j,1,n)-phi00(ne,n)) & +aimag(potential(j,1,n)-phi00(ne,n)) & *aimag(potential(j,1,n)-phi00(ne,n)) do m=2,md czn(ne,j)=czn(ne,j)+0.5*( & real(potential(j,m,n))*real(potential(j0,m,n)) & +aimag(potential(j,m,n))*aimag(potential(j0,m,n))) cznn(ne,j)=cznn(ne,j)+0.5*( & real(potential(j,m,n))*real(potential(j,m,n)) & +aimag(potential(j,m,n))*aimag(potential(j,m,n))) enddo cz(ne,j)=cz(ne,j)+real(potential(j,1,n)) czc(ne,j)=czc(ne,j) & +real(potential(j,1,n))*real(potential(j0,1,n)) & +aimag(potential(j,1,n))*aimag(potential(j0,1,n)) czcn(ne,j)=czcn(ne,j) & +real(potential(j,1,n))*real(potential(j,1,n)) & +aimag(potential(j,1,n))*aimag(potential(j,1,n)) do m=2,md cz(ne,j)=cz(ne,j)+real(potential(j,m,n)) czc(ne,j)=czc(ne,j)+0.5*( & real(potential(j,m,n))*real(potential(j0,m,n)) & +aimag(potential(j,m,n))*aimag(potential(j0,m,n))) czcn(ne,j)=czcn(ne,j)+0.5*( & real(potential(j,m,n))*real(potential(j,m,n)) & +aimag(potential(j,m,n))*aimag(potential(j,m,n))) enddo enddo enddo c c fully incorporate periodicity conditions in extensions, c where some modes have nothing to connect to. c if (iperiod.eq.2) then n_max=(nd-1)/2 do n=1,nd do nr=-n_max,n_max if (nrr(n).eq.nr) ninv(nr)=n enddo enddo do j=1,l_left-1 do m=2,md mr=mrr(m,1) beta1=2.*mr*xp*n0*qsf nshift=2.*mr*n_max*(xp+.0001)/z0 j1=j+l_right-l_left do n=-n_max-nshift,-n_max-1 if (n+nshift.le.n_max) then n1=ninv(n+nshift) bcpot=potential(j1,m,n1)*exp(-ii*beta1) cz(ne,j)=cz(ne,j)+real(bcpot) czcn(ne,j)=czcn(ne,j)+0.5*( & real(bcpot)*real(bcpot) & +aimag(bcpot)*aimag(bcpot)) endif enddo enddo enddo do j=l_right,ld do m=2,md mr=mrr(m,1) beta1=2.*mr*xp*n0*qsf nshift=2.*mr*n_max*(xp+.0001)/z0 j1=j-l_right+l_left do n=n_max+1,n_max+nshift if (n-nshift.ge.-n_max) then n1=ninv(n-nshift) bcpot=potential(j1,m,n1)*exp(ii*beta1) cz(ne,j)=cz(ne,j)+real(bcpot) czcn(ne,j)=czcn(ne,j)+0.5*( & real(bcpot)*real(bcpot) & +aimag(bcpot)*aimag(bcpot)) endif enddo enddo enddo endif c time correlation function, restarted every ninterv nprnts if (mod(ne-1,ninterv).eq.0) then l=ldb/2+1 do m=1,md do n=1,nd phit0(m,n)=potential(l,m,n) enddo enddo endif ct(ne)=0. j=ldb/2+1 do n=1,nd ct(ne)=ct(ne) & +real(potential(j,1,n))*real(phit0(1,n)) & +aimag(potential(j,1,n))*aimag(phit0(1,n)) do m=2,md ct(ne)=ct(ne)+0.5*( & real(potential(j,m,n))*real(phit0(m,n)) & +aimag(potential(j,m,n))*aimag(phit0(m,n))) enddo enddo init=0 return end