module diagnostics ! ! (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 contains subroutine prstep(t_start) ! ! Prints out various quantities of interest. ! Checks energy conservation, only valid for eps=0. ! This calculates the analytic formula ! for dW/dt, and compares with W(t+dt)-W(t-dt). ! use constants, only: ii, pi use itg_data use fields use gryffin_grid use linear use gryffin_layouts use mp, only: broadcast, proc0 use diagnostics_arrays real, dimension(m_low:m_alloc, & n_low:n_alloc) :: dummy, dum1k, & dum2k, den2m, v2m, tpar2m, tp2m, pr2m, d1eta2, d2eta2, d3eta2, & d1par2, d2par2, wupsi2, d1r2, d3r2, d4r2, d5r2, d6r2, d7r2, d8r2, & d9r2, fluxaemn ! (m,n) magnetic flutter contribution to ion particle flux real, dimension (m_low:m_alloc, & n_low:n_alloc, nspecies) :: fluxaimn ! total magnetic flutter contribution to ion particle flux real, dimension(nspecies) :: fluxai real, dimension(l_left:l_right) :: dz complex, dimension(m_low:m_alloc, n_low:n_alloc) :: phit0 real phirms1, omegad real dz0,znorm, t_start, deltat real dum,nuur,nuui,den2,v2,tpar2,tp2,vol real d1eta,d2eta,d3eta,d1par,d2par,wupsi real d1r,d3r,d4r,d5r,d6r,d7r,d8r,d9r complex, dimension(ld, m_low:m_alloc, & n_low:n_alloc, nspecies) :: totpr, & ky_psiv, nbar, ky_psiq, p_perp, kyn1, kyt1, kyt2, kwn1, kwt1, kwt2, & kpt1, kpt2, ktu1, ktt1, ktt2, ktr1, ktt3, ktt4, & ky_aparv, ubar complex, dimension(ld, m_low:m_alloc, & n_low:n_alloc) :: totpr_e, u_eeff, ky_phi, ky_apar complex, dimension(kdpass, m_low:m_alloc, & n_low:n_alloc) :: ky_phiba complex bcpot integer, dimension(-nd:nd) :: ninv integer :: init = 1 integer :: l_right_loc integer nshift, n_max, n1, j1 integer i, l, m, n, j, j0, mr, nr, mpts, npts real ky, flux1, flux2, x, y, beta1, bz, asq1, nesq1, uesq1 real fluxae logical skip_correlations totpr = 0.; ky_psiv = 0.; nbar = 0.; ky_psiq = 0.; p_perp = 0. kyn1 = 0.; kyt1 = 0.; kyt2 = 0.; kwn1 = 0.; kwt1 = 0.; kwt2 = 0. kpt1 = 0.; kpt2 = 0.; ktu1 = 0.; ktt1 = 0.; ktt2 = 0.; ktr1 = 0. ktt3 = 0.; ktt4 = 0.; ky_phi = 0. if(allocated(fluxe)) fluxe(ne) = 0. if(allocated(fluxemn)) fluxemn(ne,:,:) = 0. if(allocated(qfluxe)) qfluxe(ne) = 0. if(allocated(qfluxemn)) qfluxemn(ne,:,:) = 0. if (nparmom == 3) then nuur=nu5r nuui=nu5i else nuur=0.0 nuui=0.0 endif if (idx_local(m_lo, md, nd)) dtold = dt(md) call broadcast(dtold, proc_id(m_lo, md, nd)) ! new growth rate stuff , mbeer if(init == 1) phisqold = 0. if(init == 0) phisqold(:,:) = phisq(ne-1,:,:) call volsq(potential,potential,phirms1,dummy) do m = m_low, m_high do n = n_low, n_high phisq(ne,m,n) = dummy(m,n) enddo enddo phirms(ne) = sqrt(phirms1) do n=n_low, n_high do m=m_low, m_high if ((phisqold(m,n) /= 0.) .and. (init == 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 ! ! end of growth rate stuff ! ! Evaluate drive terms again to make dE/dt 2nd order in dt ! ! ! Driving terms ! ! ********************* ! ! load arrays to be volume averaged ! call flr(potential) do i=1,nspecies do n=n_low, n_high do m=m_low, m_high do l=1,ld ! Implement omega_gb and omega_kap correctly here at a later time. ! For now, simply use half the sum (reducing to the old usage). ! omegad=omega_gb(l,m,n,i)+omega_kap(l,m,n,i) ! 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) ! toroidal drive kwn1(l,m,n,i)=-ii*omegad*phi_nd(l,m,n,i) kwt1(l,m,n,i)=-2.*ii*omegad* & (nu2i*t_perp(l,m,n,i)+0.5*phi_tpard(l,m,n,i) ) kwt2(l,m,n,i)=-2.*ii*omegad* & ( nu3i*t_par(l,m,n,i)+0.5*phi_tperpd(l,m,n,i)) ! 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) ! toroidal damping ktu1(l,m,n,i)=2.*nu5r*abs(omegad)*u_par(l,m,n,i) ktt1(l,m,n,i)=2.*nu1r*abs(omegad)*t_par(l,m,n,i) ktt2(l,m,n,i)=2.*nu2r*abs(omegad)*t_perp(l,m,n,i) ktt3(l,m,n,i)=2.*nu3r*abs(omegad)*t_par(l,m,n,i) ktt4(l,m,n,i)=2.*nu4r*abs(omegad)*t_perp(l,m,n,i) ! 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, 0.) call kparfft(kpt2, 1, 0.) call kparfft(ktr1,-1, 0.) ! ! Do the volume averages, add to drive ! ! BD igeo=1 changes, 6/25/95: ! I haven't worked through this section carefully. It may be that I ! am wrong to replace volsq calls with volflux calls here. Since ! this is just a diagnostic, I will come back to it later. ! ! olddrive=drive ! call volflux(density(:,:,:,1),kyn1(:,:,:,1),d1eta,d1eta2) ! call volflux(t_par(:,:,:,1),kyt2(:,:,:,1),d2eta,d2eta2) ! call volflux(t_perp(:,:,:,1),kyt1(:,:,:,1),d3eta,d3eta2) ! ! ?? BD 6/25/95: ! Not generalized for igeo=1 call volsq(u_par(:,:,:,1),ktr1(:,:,:,1),wupsi,wupsi2) ! ! ?? BD 6/25/95: ! Not generalized for igeo=1 call volsq(t_par(:,:,:,1),kpt1(:,:,:,1),d1par,d1par2) ! ! ?? BD 6/25/95: ! Not generalized for igeo=1 call volsq(t_perp(:,:,:,1),kpt2(:,:,:,1),d2par,d2par2) ! call volflux(density(:,:,:,1),kwn1(:,:,:,1),d1r,d1r2) ! call volflux(t_par(:,:,:,1),kwt1(:,:,:,1),d4r,d4r2) ! call volflux(t_perp(:,:,:,1),kwt2(:,:,:,1),d7r,d7r2) ! call volflux(u_par(:,:,:,1),ktu1(:,:,:,1),d3r,d3r2) ! nu**r call volflux(t_par(:,:,:,1),ktt1(:,:,:,1),d5r,d5r2) ! nu**r call volflux(t_par(:,:,:,1),ktt2(:,:,:,1),d6r,d6r2) ! nu**r call volflux(t_perp(:,:,:,1),ktt4(:,:,:,1),d9r,d9r2) ! nu**r call volflux(t_perp(:,:,:,1),ktt3(:,:,:,1),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 ! ! ******************* ! Check energy conservation for 3+1 ! ! call volsq(density(:,:,:,1),density(:,:,:,1),den2,den2m) call volsq(u_par(:,:,:,1),u_par(:,:,:,1),v2,v2m) call volsq(t_par(:,:,:,1),t_par(:,:,:,1),tpar2,tpar2m) call volsq(t_perp(:,:,:,1),t_perp(:,:,:,1),tp2,tp2m) wenr=0.5*(den2+v2+tpar2/2.+tp2) wenx(ne)=wenr do n = n_low, n_high do m=m_low, m_high 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 ! ! Total growth rate ! if(wenr /= 0.) then if(idx_local(m_lo, md, nd)) gamx(ne)=(wenr-pwenr)/(dt(md)*wenr) call broadcast(gamx(ne), proc_id(m_lo, md, nd)) endif ! ! **************** ! ! This is the fastest way to code the new flux. ! !********************************************************************* ! Calculate the heat flux, using the full pressure ! where P = n (T_par/2 + T_perp). Thus, in perturbed quantities, ! P_1 = 3/2 n_1 + T_perp1 + 1/2 T_par1 ! Define Pressures ! if (iflr /= 6) then totpr = t_perp + 0.5*t_par + 1.5*density p_perp = density + t_perp do n = n_low, n_high do m = m_low, m_high ky=mrr(m,n)/y0 ky_psiv(:,m,n,:)=ii*ky*phi_u(:,m,n,:) ky_psiq(:,m,n,:)=ii*ky*phi_qperp(:,m,n,:) enddo enddo ! ! Do the volume averages ! if (beta_e > 0.) then call volflux (apar,apar,asq1,dummy) call volflux (n_e,n_e,nesq1,dummy) call volflux (u_e,u_e,uesq1,dummy) asq(ne)=asq1 nesq(ne)=nesq1 uesq(ne)=uesq1 endif do i=1,nspecies call volflux(totpr(:,:,:,i), ky_psiv(:,:,:,1), flux1, dummy) do n = n_low, n_high do m = m_low, m_high qfluximn(ne,m,n,i)=dummy(m,n)*n_I(i)*tau(i)*tiovte enddo enddo call volflux(p_perp(:,:,:,i), ky_psiq(:,:,:,1), flux2, dummy) do n = n_low, n_high do m = m_low, m_high 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(idx_local(m_lo, md, nd)) then if(ne == 1) then deltat=time(md)-t_start else deltat=time(md)-timo(ne-1) endif endif call broadcast(deltat, proc_id(m_lo, md, nd)) if (ifilter > 0) then if (idx_local(m_lo, 1, 1)) then chi_int = chi_int - wpfx(ne,1)*deltat/eta(1)+1.5*rdiff(1,1)*deltat endif call broadcast(chi_int, proc_id(m_lo, 1, 1)) endif else ! ! Ron's way ! do i=1,nspecies do n=n_low, n_high do m=m_low, m_high ky=mrr(m,n)/y0 totpr(:,m,n,i)=t_perp(:,m,n,i) + 0.5*t_par(:,m,n,i) + 1.5*density(:,m,n,i) & -1.5*potential(:,m,n) + phi_n(:,m,n,i)/etai & +(1.5-1./etai)*phi_u(:,m,n,i) ky_phi(:,m,n)=ii*ky*potential(:,m,n) enddo enddo enddo ! ! Do the volume averages ! call volflux(totpr(:,:,:,1),ky_phi,flux1,dummy) wpfx(ne,1)=flux1 endif ! ! Particle and electron heat fluxes ! ky_apar = 0. ky_aparv = 0. ubar = 0. totpr_e = 0. if (beta_e > 0.) then do i=1,nspecies do n=n_low, n_high do m=m_low, m_high ky=mrr(m,n)/y0 totpr_e(:,m,n)=tperp_e(:,m,n) + 0.5*tpar_e(:,m,n)+1.5*n_e(:,m,n) ubar(:,m,n,i)= (u_par(:,m,n,i)*nwgt(:,m,n,i) & +q_perp(:,m,n,i)*twgt(:,m,n,i))*g0wgtp(:,m,n,i) ky_aparv(:,m,n,i)=ii*ky*J0apar(:,m,n,i) ky_apar(:,m,n)=ii*ky*apar(:,m,n) enddo enddo enddo endif do i=1,nspecies do n=n_low, n_high do m=m_low, m_high ky=mrr(m,n)/y0 nbar(:,m,n,i)=((density(:,m,n,i)*nwgt(:,m,n,i) & +t_perp(:,m,n,i)*twgt(:,m,n,i))*g0wgtp(:,m,n,i) & -pfilter3(:,m,n,i)*potential(:,m,n))/charge(i) ky_phi(:,m,n)=ii*ky*potential(:,m,n) enddo enddo enddo do i=1,nspecies call volflux(nbar(:,:,:,i),ky_phi,flux1,dummy) fluxi(ne,i)=flux1 do n=n_low, n_high do m = m_low, m_high fluximn(ne,m,n,i) = dummy(m,n) enddo enddo enddo do n=n_low, n_high do m=m_low, m_high if(m == 1) cycle 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 if(idx_local(m_lo, m, n)) then ky=mrr(m,n)/y0 do l=1,kdpass ky_phiba(l,m,n)=ii*ky*phi_ba(l,m,n) enddo endif enddo ! ?? BD 6/25/95: ! Not generalized for igeo=1 (kappa integrals ignored so far in changes) if(epse > 0) then call volkflux(e_density,ky_phiba,flux1,dum1k) call volkflux(phi_ba,ky_phiba,flux2,dum2k) fluxe(ne)=flux1-flux2*tiovte do n=n_low, n_high do m = m_low, m_high fluxemn(ne,m,n) = dum1k(m,n) - tiovte*dum2k(m,n) enddo enddo ! ?? BD 6/25/95: ! Not generalized for igeo=1 call volkflux(e_p,ky_phiba,flux1,dum1k) qfluxe(ne)=1.5*(flux1-flux2*tiovte) do n = n_low, n_high do m = m_low, m_high qfluxemn(ne,m,n) = 1.5*(dum1k(m,n) - tiovte*dum2k(m,n)) enddo enddo endif ! PS 12/97 first stab at passing electron flux for finite beta u_eeff = 0. if (beta_e > 0.) then call volflux(n_e,ky_phi,flux1,dummy) fluxe(ne)=fluxe(ne)+flux1 ! PS!!! note using same array for em electron flux and bounce averaged ! electron flux. may want to change for finite beta w/trapped e's do n = n_low, n_high do m = m_low, m_high fluxemn(ne,m,n)=fluxemn(ne,m,n)+dummy(m,n) enddo enddo do n = n_low, n_high do m = m_low, m_high ky=mrr(m,n)/y0 u_eeff(:,m,n)=ii*ky*tpar_e(:,m,n)*n_e(:,m,n) enddo enddo call volflux(n_e, u_eeff, flux1, dummy) if(proc0) write(*,*) 'ne ,ky ne tpar_e= ',flux1 endif ! PS 12/97 first stab at passing electron heat flux for finite beta if (beta_e > 0.) then call volflux(totpr_e,ky_phi,flux1,dummy) qfluxe(ne)=1.5*flux1 do n = n_low, n_high do m = m_low, m_high qfluxemn(ne,m,n)=qfluxemn(ne,m,n)+dummy(m,n) enddo enddo endif ! PS 3/98 first stab at magnetic flutter ion particle flux if (beta_e > 0.) then fluxai = 0. fluxaimn = 0. do i=1,nspecies call volflux(ubar(:,:,:,i),ky_aparv(:,:,:,i),flux1,dummy) fluxai(i)=flux1 if (proc0) write(*,*) 'method 1 ubar,ky_aparv fluxai= ',flux1 call volflux(ubar(:,:,:,i),ky_apar,flux1,dummy) if (proc0) write(*,*) 'method 2 ubar,ky_apar fluxai= ',flux1 call volflux(u_par(:,:,:,i),ky_aparv(:,:,:,i),flux1,dummy) if (proc0) write(*,*) 'method 3 u_par,ky_aparv fluxai= ',flux1 call volflux(u_par(:,:,:,i),ky_apar,flux1,dummy) if (proc0) write(*,*) 'method 4 u_par,ky_apar fluxai= ',flux1 do n = n_low, n_high do m = m_low, m_high fluxaimn(m,n,i)=dummy(m,n) enddo enddo enddo endif ! PS 3/98 first stab at magnetic flutter electron particle flux if (beta_e > 0.) then fluxae = 0. fluxaemn = 0. call volflux(u_e, ky_apar,flux1,dummy) fluxae=flux1 if (proc0) write(*,*) 'Method 1 u_e,ky_apar fluxae= ',flux1 ! do n = n_low, n_high ! do m = m_low, m_high ! do l=1,ld-1 ! u_eeff(l,m,n)=ubar(l,m,n,1)-rkperp2(l,m,n)*apar(l,m,n)/(tiovte*beta_e) ! enddo ! enddo ! enddo ! call volflux(u_eeff, ky_apar,flux1,dummy) ! if (proc0) write(*,*) 'Method 2 u_eeff,ky_apar fluxae= ',flux1 do n = n_low, n_high do m = m_low, m_high fluxaemn(m,n)=dummy(m,n) enddo enddo call volflux(apar,ky_apar,flux1,dummy) if (proc0) write(*,*) 'apar, kyapar= ',flux1 do n = n_low, n_high do m = m_low, m_high u_eeff(:,m,n)=-rkperp2(:,m,n)*apar(:,m,n)/(tiovte*beta_e) enddo enddo call volflux(u_eeff, ky_apar,flux1,dummy) if (proc0) write(*,*) '-kperp^2 apar/(tiovte beta_e) ,ky_apar= ',flux1 endif ! temporary output PS 3/98 if (beta_e > 0. .and. proc0) write(*,2000) ne,wpfx(ne,1),qfluxe(ne), & fluxe(ne),fluxi(ne,1),fluxae,fluxai(1) 2000 format(' ne: ',i5,' Q_i: ',1pe10.2,' Q_e: ',1pe10.2, & ' Gamma_e: ',1pe10.2,' Gamma_i: ',1pe10.2,' G_flutter_e: ',1pe10.2, & ' G_flutter_i: ',1pe10.2) ! ******************* ! ! total energy, Kinetic energy, and pressure vs. mode ! (Not really kinetic energy, just v^2 ) ! ! Define total pressure ! call volsq(totpr(:,:,:,1),totpr(:,:,:,1),dum,pr2m) do n = n_low, n_high do m = m_low, m_high 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 ! ! Phi* Phi for each m and n at z=0 ! l=ldb/2 do n = n_low, n_high do m = m_low, m_high psp(ne,m,n)=real(potential(l,m,n))**2+aimag(potential(l,m,n))**2 enddo enddo ! ! Other variables we want to save ! if(idx_local(m_lo, md, nd)) timo(ne)=time(md) call broadcast(timo(ne), proc_id(m_lo, md, nd)) ! ! Flux surface averaged poloidal and toroidal flows, for m=0,n=0 mode ! if(iperiod /= 2) then l_left=1 l_right_loc=ld-1 else l_right_loc=l_right endif dz0=2.*x0/ldb do l=l_left,l_right_loc ! backwards compatibility issue: if(igeo == 0) then dz(l)=dz0*jacobian(l) else ! BD: ! factor of jacobian missing for igeo=0, which I am including ! to improve the HRW work. ! 4.6.98 dz(l)=dz0*jacobian(l) endif enddo utor(ne)=0.0 upol(ne)=0.0 phih0(ne)=0.0 denh0(ne)=0.0 uparh0(ne)=0.0 tparh0(ne)=0.0 tperph0(ne)=0.0 qparh0(ne)=0.0 qperph0(ne)=0.0 uparph0(ne)=0.0 veh0(ne)=0.0 ! ! ?? BD 6/25/95: ! Not generalized for igeo=1 ! (bz not upgraded for general geometry yet, so be careful) ! bz=1/sqrt(1+eps_over_q**2) if(idx_local(m_lo, 1, 1)) then do i=1,nspecies znorm=0. do l=l_left,l_right_loc utor(ne)=utor(ne)+dz(l)*(real(u_par(l,1,1,i)) & +eps_over_q*bmaginv(l)*bz*rkx(l,1,1)*aimag(phi_u(l,1,1,i))) & *bmaginv(l)*bz upol(ne)=upol(ne)+dz(l)*(real(u_par(l,1,1,i))*eps_over_q & -bmaginv(l)*bz*rkx(l,1,1)*aimag(phi_u(l,1,1,i))) & *bmaginv(l)*bz phih0(ne)=phih0(ne)+dz(l)*phi_u(l,1,1,i) denh0(ne)=denh0(ne)+dz(l)*density(l,1,1,i) uparh0(ne)=uparh0(ne)+dz(l)*u_par(l,1,1,i) tparh0(ne)=tparh0(ne)+dz(l)*t_par(l,1,1,i) tperph0(ne)=tperph0(ne)+dz(l)*t_perp(l,1,1,i) qparh0(ne)=qparh0(ne)+dz(l)*q_par(l,1,1,i) qperph0(ne)=qperph0(ne)+dz(l)*q_perp(l,1,1,i) uparph0(ne)=uparph0(ne)+dz(l)*u_par(l,1,1,i)*eps_over_q veh0(ne)=veh0(ne)+dz(l)*phi_u(l,1,1,i)*ii*rkx(l,1,1)/b_mag(l) znorm=znorm+dz(l) enddo enddo if(znorm == 0.) write(*,*) 'znorm = 0: problem.' utor(ne)=utor(ne)/znorm upol(ne)=upol(ne)/znorm phih0(ne)=phih0(ne)/znorm denh0(ne)=denh0(ne)/znorm uparh0(ne)=uparh0(ne)/znorm tparh0(ne)=tparh0(ne)/znorm tperph0(ne)=tperph0(ne)/znorm qparh0(ne)=qparh0(ne)/znorm qperph0(ne)=qperph0(ne)/znorm uparph0(ne)=uparph0(ne)/znorm veh0(ne)=veh0(ne)/znorm endif call broadcast(utor, proc_id(m_lo, 1, 1)) call broadcast(upol, proc_id(m_lo, 1, 1)) call broadcast(phih0, proc_id(m_lo, 1, 1)) call broadcast(denh0, proc_id(m_lo, 1, 1)) call broadcast(uparh0, proc_id(m_lo, 1, 1)) call broadcast(tparh0, proc_id(m_lo, 1, 1)) call broadcast(tperph0, proc_id(m_lo, 1, 1)) call broadcast(qparh0, proc_id(m_lo, 1, 1)) call broadcast(qperph0, proc_id(m_lo, 1, 1)) call broadcast(uparh0, proc_id(m_lo, 1, 1)) call broadcast(veh0, proc_id(m_lo, 1, 1)) ! write(*,*) 've/ve0: ',veh0(ne)/veh0(1) ! ! flux surface averaged quantities phi00(ne,:)=0.0 den00(ne,:)=0.0 upar00(ne,:)=0.0 tpar00(ne,:)=0.0 tperp00(ne,:)=0.0 do n=n_low, n_high m = 1 if(idx_local(m_lo, m, n)) then znorm=0. do l=l_left,l_right_loc-1 phi00(ne,n)=phi00(ne,n)+potential(l,m,n)*dz(l) den00(ne,n)=den00(ne,n)+density(l,m,n,1)*dz(l) upar00(ne,n)=upar00(ne,n)+u_par(l,m,n,1)*dz(l) tpar00(ne,n)=tpar00(ne,n)+t_par(l,m,n,1)*dz(l) tperp00(ne,n)=tperp00(ne,n)+t_perp(l,m,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 endif enddo n_e00(ne,:)=0.0 u_e00(ne,:)=0.0 apar00(ne,:)=0.0 if (beta_e > 0.) then do n=n_low, n_high m = 1 if(idx_local(m_lo, m, n)) then znorm=0. do l=l_left,l_right_loc-1 n_e00(ne,n)=n_e00(ne,n)+n_e(l,m,n)*dz(l) u_e00(ne,n)=u_e00(ne,n)+u_e(l,m,n)*dz(l) apar00(ne,n)=apar00(ne,n)+apar(l,m,n)*dz(l) znorm=znorm+dz(l) enddo n_e00(ne,n)=n_e00(ne,n)/znorm u_e00(ne,n)=u_e00(ne,n)/znorm apar00(ne,n)=apar00(ne,n)/znorm endif enddo endif ! Only radial correlation functions in parallel code, at least for now. ! BD/QL 10/21/98 ! Correlation functions skip_correlations=.true. if(skip_correlations) goto 20 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=m_low, m_high ! not correct if kx's are distributed over processors: ?? do n=n_low, n_high nr=nrr(n) cx(ne,j)=cx(ne,j)+real(potential(l,m,n)*exp(ii*nr*x)) enddo enddo enddo mpts=4*md 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 mr=mrr_all(m,n) cy(ne,j)=cy(ne,j)+real(potential(l,m,n)*exp(-ii*mr*y)) enddo enddo enddo ! ??BD 9/22/95: ! Apparently these integrals need Jacobian factors if igeo=1, ! but I am leaving it for later: 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 ! ! fully incorporate periodicity conditions in extensions, ! where some modes have nothing to connect to. ! if (iperiod == 2) then n_max=(nd-1)/2 do n=1,nd do nr=-n_max,n_max if (nrr(n) == nr) ninv(nr)=n enddo enddo do j=1,l_left-1 do m=2,md mr=mrr_all(m,1) beta1=2.*mr*xp*n0*qsf nshift=2.*mr*n_max*(xp+.0001)/z0 j1=j+l_right_loc-l_left do n=-n_max-nshift,-n_max-1 if (n+nshift <= 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_loc,ld do m=2,md mr=mrr_all(m,1) beta1=2.*mr*xp*n0*qsf nshift=2.*mr*n_max*(xp+.0001)/z0 j1=j-l_right_loc+l_left do n=n_max+1,n_max+nshift if (n-nshift >= -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 ! time correlation function, restarted every ninterv nprnts phit0 = 0. if (mod(ne-1,ninterv) == 0) then l=ldb/2+1 phit0(:,:)=potential(l,:,:) 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 20 continue init=0 end subroutine prstep subroutine preenr ! ! (* old: Determine wkin,wmgt,wki and wenr at time step=ne-1 *) ! ! ! Calculates energy (wner) at time step=ne-1 ! Also stores grtmx (mode amplitudes ! use fields, only: density, u_par, t_par, t_perp, potential use itg_data, only: ne use gryffin_layouts use diagnostics_arrays real, dimension(m_low:m_alloc, & n_low:n_alloc) :: dummy, u2 real dum integer nea, m, n real den2,v2,tpar2,tp2 ! call volsq(density(:,:,:,1),density(:,:,:,1),den2,dummy) call volsq(u_par(:,:,:,1),u_par(:,:,:,1),v2,dummy) call volsq(t_par(:,:,:,1),t_par(:,:,:,1),tpar2,dummy) call volsq(t_perp(:,:,:,1),t_perp(:,:,:,1),tp2,dummy) call volsq(potential,potential,dum,u2) ! if (nparmom == 3) then ! nuur=nu5r ! nuui=nu5i ! else ! nuur=0.0 ! nuui=0.0 ! endif nea=ne+1 do n = n_low, n_high do m = m_low, m_high grtmx(nea,m,n) = u2(m,n) enddo enddo pwenr=0.5*(den2+v2+tpar2/2.+tp2) end subroutine preenr subroutine posten ! ! (* old comment Determine wkin,wmgt,wki and wenr at time step=ne-1 ! like preenr but had an option to adjust dt. *) ! ! ! I'm not sure why, it probably has something to do with variable ! timestep, but here grtmx is recalculated. ! ! use mp, only: broadcast use fields, only: potential, density, u_par, t_par, t_perp use itg_data, only: ncycle, dt, ne use gryffin_layouts use gryffin_grid, only: md, nd use diagnostics_arrays real, dimension(m_low:m_alloc, & n_low:n_alloc) :: u2, dummy real edif,etot,dum, dt1 real den2,v2,tpar2,tp2,wenrn,dedt integer m,n call volsq(potential,potential,dum,u2) if(ncycle > 2) then do n=n_low, n_high do m=m_low, m_high edif=u2(m,n)-grtmx(ne,m,n) etot=u2(m,n)+grtmx(ne,m,n) if (abs(etot) < 1.0e-33) then grtmx(ne,m,n) = 0.0 else grtmx(ne,m,n)=edif/(2.*dt(m)*etot) endif enddo enddo else do n=n_low, n_high do m = m_low, m_high grtmx(ne,m,n)=0.0 enddo enddo endif ! ! calculate energy conservation, 2nd order accurate in dt ! call volsq(density(:,:,:,1),density(:,:,:,1),den2,dummy) call volsq(u_par(:,:,:,1),u_par(:,:,:,1),v2,dummy) call volsq(t_par(:,:,:,1),t_par(:,:,:,1),tpar2,dummy) call volsq(t_perp(:,:,:,1),t_perp(:,:,:,1),tp2,dummy) ! ! Without subcycling, this calculation only makes sense if all modes are ! being advanced with the same time step. It is not vital in the ! linear cases, so replace dt with dt(md) as elsewhere: ! if (idx_local(m_lo, md, nd)) dt1 = dt(md) call broadcast(dt1, proc_id(m_lo, md, nd)) wenrn=0.5*(den2+v2+tpar2/2.+tp2) dedt=(dt1/dtold*(wenr-pwenr)+dtold/dt1*(wenrn-wenr))/(dt1+dtold) if(wenrn /= 0.) then pcerr(ne)=(dedt-drive)/wenrn else pcerr(ne)=0. endif end subroutine posten subroutine volsq(a,b,axb,axb_by_mode) ! ! Volume integral of q*p / Volume ! ! wp is volume integral over ! (-zlen/2,zlen/2)(-ylen/2,ylen/2)(-xlen/2,xlen/2) ! wpp is z integral for each m,n ! ! where zlen is 2.*xp if iperiod=2 and 2.*x0 otherwise. ! ! BD 6/22/95: ! dVol picks up a factor of jacobian in general geometry. ! This would change our previous answers slightly, since we ! ignored 1/B(theta) in the previous volume averages. To recover ! the earlier results, the variable "jacobian" should be independent ! of theta for igeo=0. However, this doesn't work because the ! jacobian factor in poisson retained its theta dependence. Thus, ! the logic here may look a little odd: ignore the jacobian factor ! if igeo=0, (to compare with old results) but not for igeo=1, even ! if itor=0. ! ! This means that igeo=1, itor=0 nonlinear runs will not give same ! result as igeo=0 runs! ! use mp, only: sum_allreduce use itg_data, only: igeo, jacobian, iperiod use gryffin_layouts use gryffin_grid, only: ld, ldb, md, l_left, l_right, x0 use diagnostics_arrays complex, dimension(:,m_low:,n_low:) :: a, b real, dimension(m_low:,n_low:) :: axb_by_mode real, dimension(l_left:l_right) :: dz real, dimension(md) :: zero ! integer :: init = 1 real axb,vnorm,dz0 integer l,m,n,l_right_loc dz0=2.*x0/ldb if(iperiod /= 2) then l_left=1 l_right_loc=ld-1 else l_right_loc=l_right endif ! ylen=2.*pi*y0 ! xlen=2.*pi*y0*(float(nd-1)/2.)/(shr*z0) ! ! This logic is wrong for input md=-1. ! zero(1)=1.0 do m=2,md zero(m)=0.5 enddo do l=l_left,l_right_loc if(igeo == 0) then dz(l)=dz0 else dz(l)=dz0*jacobian(l) endif enddo axb=0.0 axb_by_mode = 0. do n = n_low, n_high do m = m_low, m_high vnorm=0. do l=l_left,l_right_loc vnorm=vnorm+dz(l) axb_by_mode(m,n)=axb_by_mode(m,n) & +(real(a(l,m,n))*real(b(l,m,n)) & +aimag(a(l,m,n))*aimag(b(l,m,n)))*dz(l) enddo axb_by_mode(m,n)=axb_by_mode(m,n)/vnorm*zero(m) axb=axb+axb_by_mode(m,n) enddo enddo call sum_allreduce(axb) ! ! Debug: ! if(init == 1) write(*,*) 'int(J dtheta)= ',vnorm ! if(init == 1) init=0 ! end subroutine volsq subroutine volsqk(a,b,axb,axb_by_mode) ! Volume integral of q*p, where a and b are pitch angle moments ! ! wp is volume integral over ! (-zlen/2,zlen/2)(-ylen/2,ylen/2)(-xlen/2,xlen/2) ! wpp is z integral for each m,n ! ! ??BD 9/22/95: ! Not yet generalized for igeo=1. (I haven't checked the kappa dependence.) use mp, only: sum_allreduce use itg_data, only: iperiod, igeo, jacobian use gryffin_layouts use gryffin_grid, only: ld, ldb, md, l_left, l_right, kd, kdpass, x0 use bounce_avg, only: dkap, wgtba2 use diagnostics_arrays complex, dimension(:,m_low:,n_low:) :: a, b real, dimension(m_low:,n_low:) :: axb_by_mode real, dimension(l_left:l_right) :: dz real, dimension(md) :: zero real, dimension(kdpass, m_low:m_alloc, & n_low:n_alloc) :: c, d real axb, vnorm integer l,m,n,j,l0,l_right_loc if(iperiod /= 2) then l_left = 1 l_right_loc=ld-1 else l_right_loc=l_right endif do l=l_left,l_right_loc if(igeo == 0) then dz(l)=2.*x0/ldb else dz(l)=2.*x0/ldb*jacobian(l) endif enddo zero(1)=1.0 do m=2,md ! zero(m)=2.0 zero(m)=0.5 enddo ! ! first multiply pitch angle functions together ! do n=n_low, n_high do m=m_low, m_high if (m == 1) cycle do l=1,kd c(l,m,n)=real(a(l,m,n))*real(b(l,m,n))+aimag(a(l,m,n))*aimag(b(l,m,n)) enddo enddo m=1 if(idx_local(m_lo, m, n)) then do l=1,kdpass c(l,m,n)=real(a(l,m,n))*real(b(l,m,n))+aimag(a(l,m,n))*aimag(b(l,m,n)) enddo endif enddo ! ! now average in kappa ! do m = m_low, m_high if (m == 1) cycle do n = n_low, n_high do j=1,ldb l0=1+abs(j-ldb/2-1) d(j,m,n)=0. do l=l0,kd d(j,m,n)=d(j,m,n)+dkap(l)*c(l,m,n)*wgtba2(l,j) enddo enddo enddo enddo m=1 do n=n_low, n_high if(idx_local(m_lo, m, n)) then do j=1,ldb l0=1+abs(j-ldb/2-1) d(j,m,n)=0. do l=l0,kd-1 d(j,m,n)=d(j,m,n)+dkap(l)*(c(l,m,n)+c(l+1,m,n))/2*wgtba2(l,j) enddo do l=max(kd,l0),kd+1 d(j,m,n)=d(j,m,n)+dkap(l)*c(l,m,n)*wgtba2(l,j) enddo do l=kd+2,kdpass d(j,m,n)=d(j,m,n)+dkap(l)*(c(l,m,n)+c(l-1,m,n))/2*wgtba2(l,j) enddo enddo endif enddo ! ! finally volume average ! axb=0.0 axb_by_mode=0. do n = n_low, n_high do m = m_low, m_high vnorm=0. do l=l_left,l_right_loc axb_by_mode(m,n)=axb_by_mode(m,n)+d(l,m,n)*dz(l)*zero(m) vnorm=vnorm+dz(l) enddo axb_by_mode(m,n)=axb_by_mode(m,n)/vnorm axb=axb+axb_by_mode(m,n) enddo enddo call sum_allreduce(axb) end subroutine volsqk subroutine volflux(a, b, axb, axb_by_mode) ! ! Volume integral of a*b*grad(rho)/volume integral of grad(rho) ! ! For a vector flux F, the quantity that is needed for transport ! calculations is dV/drho* < F.grad(rho) >, where <...> signifies a ! flux surface average: ! ! < G > == Int(dtheta dalpha drho J G) / Int(J dtheta dalpha drho) ! ! / ! == Int(dtheta dalpha drho J G) / dV/drho * delta_rho ! / ! ! Taking delta_rho -> 0 is the definition of the flux surface average. ! Our simulation domain is small enough that averaging over the entire ! radial domain is ok. ! ! Thus, dV/drho = 2.*pi*Int(J dtheta) in our coordinates, since J is ! only a function of theta. ! ! This routine actually calculates ! ! < F.grad rho > / < |grad rho| > = ! ! Int(dtheta dalpha drho J |grad rho| F_rho) / Int(J dtheta dalpha drho) ! ----------------------------------------------------------------------- ! Int(dtheta dalpha drho J |grad rho| ) / Int(J dtheta dalpha drho) ! ! F.A delta(rho) ! = -------------- ! A delta(rho) ! ! so that, for example, the heat flux that is reported is the heat ! flux per unit area. ! ! This quantity is calculated assuming that F is an ExB type flux, and ! that a and b are the parts whose product is the generalized radial ! (rho direction) component of F, denoted F_rho. ! ! Note that < |grad rho| > * dV/drho is the area of the flux surface. ! ! dV/drho == 2.*pi*Int(dtheta J(theta)), where -N*pi <= theta <= N*pi. ! dV/drho == 2.*pi*vnorm here (the 2*pi is from Int(dalpha)) ! ! This modification for igeo=1 is useful because it reduces to the ! previous definition used for the volume-averaged flux in the ! circular flux surface limit. Also, normalization issues associated ! with our box actually extending multiples of 2*pi in the theta ! direction are cleanly handled with this choice, as follows: ! ! Our reported heat flux fits into the the fundamental transport ! equation ! ! (3/2) d nT/dt + 1/V' d/drho V' < Q.grad rho> + ... = 0 ! ! like this: ! ! (3/2) d nT/ dt + 1/V' d/drho A Q_itgc + ... = 0 ! ! where Q_itgc is the heat flux calculated with this subroutine, ! and A is the area of the flux surface. In this equation, ! V' should be the physical V' (that is, using pi < theta < pi). ! ! This transport equation can be rewritten (using definitions above): ! ! (3/2) d nT/ dt + (< |grad rho| >/A) d/drho A Q_itgc + ... = 0. ! ! In this form, there is no ambiguity about which V' definition to use; ! one just calculates < |grad rho| > and the surface area, plugs in ! the heat flux calculated from itgc, and that's that. Note that ! there will be an additional factor of <|grad rho|> (typically ! stabilizing, I believe) that this definition of Q fails to capture. ! ! The GA people tend (as I understand it) to lump the factor of ! <|grad rho|> out front together with the factor in the Q definition, ! and thus talk about a <|grad rho|**2> effect on the anomalous transport ! coefficients. ! ! Note that the definitions of the gradient scale lengths as ! ! 1/L_T == - d ln(T) / d ln(rho) ! ! become functions of the choice of rho, the flux-surface label. ! If one is calculating everything within this ptor-itgc-eik code, ! there is no complication arising from this observation. One may use ! any definition available. However, otherwise one must be careful to ! keep everything consistent. Factors of drho_1/drho_2 appear... ! ! Mike K. and I suggest that since we are performing local simulations ! (in radius), the default definition of rho that we choose should be ! a locally defined quantity. Furthermore, it should be very fast to ! calculate (to keep the time spent in the eikcoefs subroutine to a ! minimum). The definition which best seems to fit these two criteria ! ! rho_d = d/D ! ! where d == the diameter of the flux surface, and D == the diameter of ! the last closed flux surface. ! ! Any interpolation formulas or databases, etc., that we develop will ! probably use this definition of rho unless someone thinks of a good ! reason to go with another choice. ! use itg_data, only: igeo, gradrho, jacobian, iperiod use gryffin_grid, only: l_left, l_right, ld, md, ldb, x0 use mp, only: sum_allreduce use gryffin_layouts use diagnostics_arrays complex, dimension(:,m_low:,n_low:) :: a, b real, dimension(m_low:,n_low:) :: axb_by_mode real, dimension(l_left:l_right) :: dz real, dimension(md) :: zero real :: axb, vnorm integer l, m, n, l_right_loc zero(1)=1.0 do m=2,md zero(m)=0.5 enddo if(iperiod /= 2) then l_left=1 l_right_loc=ld-1 else l_right_loc=l_right endif do l=l_left,l_right_loc if(igeo == 0) then dz(l)=2.*x0/ldb*gradrho(l) else dz(l)=2.*x0/ldb*jacobian(l)*gradrho(l) endif enddo axb = 0.0 axb_by_mode = 0.0 do n = n_low, n_high do m = m_low, m_high vnorm=0. do l=l_left,l_right_loc vnorm=vnorm+dz(l) axb_by_mode(m,n)=axb_by_mode(m,n) & +(real(a(l,m,n))*real(b(l,m,n)) & +aimag(a(l,m,n))*aimag(b(l,m,n)))*dz(l) enddo axb_by_mode(m,n)=axb_by_mode(m,n)/vnorm*zero(m) axb=axb+axb_by_mode(m,n) enddo enddo call sum_allreduce(axb) end subroutine volflux subroutine volkflux(a,b,axb,axb_by_mode) ! ! See comment in volflux for igeo=1 generalization. ! ! Volume integral of q*p, where a and b are pitch angle moments ! ! wp is volume integral over ! (-zlen/2,zlen/2)(-ylen/2,ylen/2)(-xlen/2,xlen/2) ! wpp is z integral for each m,n ! ! ??BD 9/22/95: ! Not yet generalized for igeo=1. (I haven't checked the kappa dependence.) use mp, only: sum_allreduce use itg_data, only: igeo, gradrho, jacobian, iperiod use bounce_avg, only: dkap, wgtba2 use gryffin_layouts use gryffin_grid, only: ld, ldb, md, l_left, l_right, kd, kdpass, x0 use diagnostics_arrays complex, dimension(:,m_low:,n_low:) :: a, b real, dimension(m_low:,n_low:) :: axb_by_mode real, dimension(l_left:l_right) :: dz real, dimension(md) :: zero real, dimension(kdpass, m_low:m_alloc, & n_low:n_alloc) :: c, d real vnorm, axb integer l,m,n,j,l0, l_right_loc if(iperiod /= 2) then l_left=1 l_right_loc=ld-1 else l_right_loc=l_right endif do l=l_left,l_right_loc if(igeo == 0) then dz(l)=2.*x0/ldb*gradrho(l) else dz(l)=2.*x0/ldb*jacobian(l)*gradrho(l) endif enddo zero(1)=1.0 do m=2,md zero(m)=0.5 enddo ! ! first multiply pitch angle functions together ! do n = n_low, n_high do m = m_low, m_high if(m == 1) cycle do l=1,kd c(l,m,n)=real(a(l,m,n))*real(b(l,m,n))+aimag(a(l,m,n))*aimag(b(l,m,n)) enddo enddo m=1 if(idx_local(m_lo, m, n))then do l=1,kdpass c(l,m,n)=real(a(l,m,n))*real(b(l,m,n))+aimag(a(l,m,n))*aimag(b(l,m,n)) enddo endif enddo ! ! now average in kappa ! do m = m_low, m_high if(m == 1) cycle do n = n_low, n_high do j=1,ldb l0=1+abs(j-ldb/2-1) d(j,m,n)=0. do l=l0,kd d(j,m,n)=d(j,m,n)+dkap(l)*c(l,m,n)*wgtba2(l,j) enddo enddo enddo enddo m=1 do n=n_low, n_high if(idx_local(m_lo, m, n)) then do j=1,ldb l0=1+abs(j-ldb/2-1) d(j,m,n)=0. do l=l0,kd-1 d(j,m,n)=d(j,m,n)+dkap(l)*(c(l,m,n)+c(l+1,m,n))/2*wgtba2(l,j) enddo do l=max(kd,l0),kd+1 d(j,m,n)=d(j,m,n)+dkap(l)*c(l,m,n)*wgtba2(l,j) enddo do l=kd+2,kdpass d(j,m,n)=d(j,m,n)+dkap(l)*(c(l,m,n)+c(l-1,m,n))/2*wgtba2(l,j) enddo enddo endif enddo ! ! finally volume average ! axb=0.0 axb_by_mode=0. do n = n_low, n_high do m = m_low, m_high vnorm=0. do l=l_left,l_right_loc axb_by_mode(m,n)=axb_by_mode(m,n)+d(l,m,n)*dz(l)*zero(m) vnorm=vnorm+dz(l) enddo axb_by_mode(m,n)=axb_by_mode(m,n)/vnorm axb=axb+axb_by_mode(m,n) enddo enddo call sum_allreduce(axb) end subroutine volkflux subroutine pcons ! ! checks particle conservation ! use itg_data, only: epse use fields, only: e_density, phi_ba, potential, e_denk use gryffin_grid, only: kdpass, md, nd, kd, x0, ldb, r use bounce_avg, only: kapav use diagnostics_arrays complex, dimension(:,:,:), allocatable :: den, dent, denp complex, dimension(:), allocatable :: ntot, nt, np real dtheta integer n,l,m allocate (den(kdpass, md, nd), dent(kdpass, md, nd), denp(kdpass, md, nd)) allocate (ntot(nd), nt(nd), np(nd)) ! ! load trapped, passing, and total density ! do n=1,nd do l=1,kd dent(l,1,n)=e_density(l,1,n)-phi_ba(l,1,n) denp(l,1,n)=0. den(l,1,n)=e_density(l,1,n)-phi_ba(l,1,n) do m=2,md dent(l,m,n)=0. denp(l,m,n)=0. den(l,m,n)=0. enddo enddo do l=kd+1,kdpass dent(l,1,n)=0. denp(l,1,n)=e_density(l,1,n)-phi_ba(l,1,n) den(l,1,n)=e_density(l,1,n)-phi_ba(l,1,n) do m=2,md dent(l,m,n)=0. denp(l,m,n)=0. den(l,m,n)=0. enddo enddo enddo call kapav(phi_ba,den) ! ! flux surface average total electron density ! dtheta=2*x0/float(ldb) do n=1,nd ntot(n)=0. do l=1,ldb ntot(n)=ntot(n)+(potential(l,1,n)+e_denk(l,1,n))*(1+epse*cos(r(l)))*dtheta enddo enddo call kapav(phi_ba,dent) ! ! flux surface average trapped electron density ! dtheta=2*x0/float(ldb) do n=1,nd nt(n)=0. do l=1,ldb nt(n)=nt(n)+(potential(l,1,n)*sqrt(epse*(1+cos(r(l)))/ & (1+epse*cos(r(l))))+e_denk(l,1,n))*(1+epse*cos(r(l)))*dtheta enddo enddo call kapav(phi_ba,denp) ! ! flux surface average passing electron density ! dtheta=2*x0/float(ldb) do n=1,nd np(n)=0. do l=1,ldb np(n)=np(n)+(potential(l,1,n)*(1.-sqrt(epse*(1+cos(r(l)))/ & (1+epse*cos(r(l)))))+e_denk(l,1,n))*(1+epse*cos(r(l)))*dtheta enddo enddo do n=1,nd write(*,101) n,ntot(n),nt(n),np(n) enddo 101 format(' n,tot,trap,pass: ',i3,3(1pe14.6,1pe14.6,'; ')) ! check that flux surface average of <1/taub>_kappa = (1+epse)/2 ! do l=1,kdpass ! denp(l,1,1)=sqrt(2*epse/(1+epse))/taub(l) ! enddo ! call kapav(phi_ba,denp) ! np(1)=0 ! do l=1,ldb ! np(1)=np(1)+e_denk(l,1,1)*(1+epse*cos(r(l)))*dtheta ! enddo ! write(*,*) '<<1/taub>_kappa>_theta, (1+epse)/2: ', ! . np(1),(1+epse)/2 deallocate (den, dent, denp, ntot, nt, np) end subroutine pcons subroutine init_diagnostics use itg_data, only: debug_plotlabel use diagnostics_arrays integer, dimension(7) :: itime character(20) :: datestamp, timestamp, zone interface subroutine x05aae(i) integer, dimension(7) :: i end subroutine x05aae end interface interface function x05abe(i) character*30 x05abe integer, dimension(7) :: i end function x05abe end interface ! ! Get the time and date of this run (used to label plots): ! datestamp(:) = ' ' timestamp(:) = ' ' zone(:) = ' ' call date_and_time (datestamp, timestamp, zone) ! convert to format "year/mm/dd hr:mm:ss +-zone" if(debug_plotlabel == ' ' ) then date = & datestamp(1:4) // "/" // datestamp(5:6) // "/" // datestamp(7:8) & // " " // timestamp(1:2) // ":" // timestamp(3:4) // ":" // & timestamp(5:6) // " " // zone else date=debug_plotlabel endif ! call x05aaE(itime) ! if(debug_plotlabel == ' ' ) then ! date=x05abE(itime) ! else ! date=debug_plotlabel ! endif end subroutine init_diagnostics end module diagnostics