module diagnostics
 
  use itg_data
  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 flr_terms
    
    real, dimension(:,:), pointer, save :: phisqold
    real, dimension(:,:), allocatable :: dummy, dum1k, dum2k, den2m, &
         v2m, tpar2m, tp2m, pr2m, d1eta2, d2eta2, d3eta2, d1par2, d2par2, &
         wupsi2, d1r2, d3r2, d4r2, d5r2, d6r2, d7r2, d8r2, d9r2
    real, dimension(:), pointer, save :: dz
    complex, dimension(:,:),   allocatable :: phit0
    
    real phirms1
    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, allocatable, dimension(:,:,:,:) :: totpr, ky_psiv, nbar, &
         ky_psiq, p_perp, kyn1, kyt1, kyt2, kwn1, kwt1, kwt2, &
         kpt1, kpt2, ktu1, ktt1, ktt2, ktr1, ktt3, ktt4
    complex, allocatable, dimension(:,:,:) :: ky_phiba
    complex bcpot
    
    integer, dimension(:), pointer, save :: ninv
    
    integer :: init = 1
    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
    
    interface 
       subroutine kparfft(a, i)
         complex, dimension(:,:,:,:) :: a
         integer :: i
       end subroutine kparfft
    end interface        

    allocate (dummy(md, nd), dum1k(md, nd), dum2k(md,nd), &
         den2m(md, nd), v2m(md, nd), tpar2m(md, nd), tp2m(md, nd), &
         pr2m(md, nd), d1eta2(md, nd), d2eta2(md, nd), d3eta2(md, nd), &
         d1par2(md, nd), d2par2(md, nd), wupsi2(md, nd), d1r2(md, nd), &
         d3r2(md, nd), d4r2(md, nd), d5r2(md, nd), d6r2(md, nd), &
         d7r2(md, nd), d8r2(md, nd), d9r2(md, nd))
    allocate (totpr(ld, md, nd, nspecies), ky_psiv(ld, md, nd, nspecies), &
         nbar(ld, md, nd, nspecies), ky_psiq(ld, md, nd, nspecies), &
         p_perp(ld, md, nd, nspecies), kyn1(ld, md, nd, nspecies), &
         kyt1(ld, md, nd, nspecies), kyt2(ld, md, nd, nspecies), &
         kwn1(ld, md, nd, nspecies), kwt1(ld, md, nd, nspecies), &
         kwt2(ld, md, nd, nspecies), kpt1(ld, md, nd, nspecies), &
         kpt2(ld, md, nd, nspecies), ktu1(ld, md, nd, nspecies), &
         ktt1(ld, md, nd, nspecies), ktt2(ld, md, nd, nspecies), &
         ktr1(ld, md, nd, nspecies), ktt3(ld, md, nd, nspecies), &
         ktt4(ld, md, nd, nspecies))
    
    if (nparmom == 3) then
       nuur=nu5r
       nuui=nu5i
    else
       nuur=0.0
       nuui=0.0
    endif
    
    dtold=dt(md)
    
!     new growth rate stuff , mbeer
  
    if(init == 1) then
       allocate (phisqold(md,nd))
       phisqold = 0.
    endif
    if(init == 0) phisqold(:,:) = phisq(ne-1,:,:)
    
    call volsq(potential,potential,phirms1,dummy)
    phisq(ne,:,:) = dummy(:,:)
    phirms(ne) = sqrt(phirms1)
    
    do n=1,nd
       do m=1,md
          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 m=1,md
          do n=1,nd
             do l=1,ld
!     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(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))
!     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(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)
!     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)
!     
!     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
! <n ..dpsi/dy>
    call volflux(density,kyn1,d1eta,d1eta2,1)
! <tpar .. dpsi/dy>
    call volflux(t_par,kyt2,d2eta,d2eta2,1)
! <tperp .. dpsi/dy>
    call volflux(t_perp,kyt1,d3eta,d3eta2,1)
! <upar dpsi/dz> 
! ?? BD 6/25/95:
! Not generalized for igeo=1
    call volsq(u_par(:,:,:,1),ktr1(:,:,:,1),wupsi,wupsi2)
! <tpar |kpar| tpar>
! ?? BD 6/25/95:
! Not generalized for igeo=1
    call volsq(t_par(:,:,:,1),kpt1(:,:,:,1),d1par,d1par2)
! <tperp |kpar| tperp>
! ?? BD 6/25/95:
! Not generalized for igeo=1
    call volsq(t_perp(:,:,:,1),kpt2(:,:,:,1),d2par,d2par2)
! <n w_d dpsi/dy>
    call volflux(density,kwn1,d1r,d1r2,1)
! <tpar w_d ... dpsi/dy>
    call volflux(t_par,kwt1,d4r,d4r2,1)
! <tperp w_d ... dpsi/dy>
    call volflux(t_perp,kwt2,d7r,d7r2,1)
! <upar |w_d k_y| upar>
    call volflux(u_par,ktu1,d3r,d3r2,1)
!  nu**r <t* |w_d k_y| t*>
    call volflux(t_par,ktt1,d5r,d5r2,1)
!  nu**r <t* |w_d k_y| t*>
    call volflux(t_par,ktt2,d6r,d6r2,1)
!  nu**r <t* |w_d k_y| t*>
    call volflux(t_perp,ktt4,d9r,d9r2,1)
!  nu**r <t* |w_d k_y| t*>
    call volflux(t_perp,ktt3,d8r,d8r2,1)

    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
  
    wenrk(ne,:,:)= 0.5*(den2m(:,:) + v2m(:,:) + tpar2m(:,:)/2. + tp2m(:,:))
    dketa(ne,:,:)= -d1eta2(:,:) - 0.5*d2eta2(:,:) - d3eta2(:,:)
    dkpar(ne,:,:)= -0.5*d1par2(:,:) - d2par2(:,:)
    dktor(ne,:,:)= -d1r2(:,:) - d4r2(:,:) - d7r2(:,:)
    dktdp(ne,:,:)= -d3r2(:,:) - d5r2(:,:)/2. - d6r2(:,:)/2. - d8r2(:,:) - d9r2(:,:)
    wkups(ne,:,:)= -wupsi2(:,:)
    dktot(ne,:,:)= dketa(ne,:,:) + dkpar(ne,:,:) + dktor(ne,:,:) + dktdp(ne,:,:) + wkups(ne,:,:)
!     
!  Total growth rate
!
    if(wenr /= 0.) gamx(ne)=(wenr-pwenr)/(dt(md)*wenr)
!
!   ****************    
!     
!     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=1,nd
          do m=1,md
             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     
!
       do i=1,nspecies
          call volflux(totpr, ky_psiv, flux1, dummy, i)
          qfluximn(ne,:,:,i)=dummy(:,:)*n_I(i)*tau(i)*tiovte

          call volflux(p_perp, ky_psiq, flux2, dummy, i)

          qfluximn(ne,:,:,i)=qfluximn(ne,:,:,i)+dummy(:,:)*n_I(i)*tau(i)*tiovte

          wpfx(ne,i)=(flux1+flux2)*n_I(i)*tau(i)*tiovte
       enddo

       if(ne == 1) then
          deltat=time(md)-t_start
       else
          deltat=time(md)-timo(ne-1)
       endif
       if(ifilter > 0) chi_int=chi_int-wpfx(ne,1)*deltat/eta(1)+1.5*rdiff(1,1)*deltat              
    else
!     
!     Ron's way
!
       do i=1,nspecies
          do n=1,nd
             do m=1,md
                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_psiv(:,m,n,i)=ii*ky*potential(:,m,n)
             enddo
          enddo
       enddo
!     
! Do the volume averages     
!
       call volflux(totpr,ky_psiv,flux1,dummy,1)

       wpfx(ne,1)=flux1
    endif
!
! Particle and electron heat fluxes
!
    do i=1,nspecies
       do n=1,nd
          do m=1,md
             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_psiv(:,m,n,i)=ii*ky*potential(:,m,n)
          enddo
       enddo
    enddo
    do i=1,nspecies
       call volflux(nbar,ky_psiv,flux1,dummy,i)
       fluxi(ne,i)=flux1
       fluximn(ne,:,:,i) = dummy(:,:)
    enddo
    
    allocate (ky_phiba(kdpass, md, nd))
    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
! ?? 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
       fluxemn(ne,:,:) = dum1k(:,:) - tiovte*dum2k(:,:)
! ?? BD 6/25/95:
! Not generalized for igeo=1
       call volkflux(e_p,ky_phiba,flux1,dum1k)
       qfluxe(ne)=1.5*(flux1-flux2*tiovte)
       qfluxemn(ne,:,:) = 1.5*(dum1k(:,:) - tiovte*dum2k(:,:))
    endif
    deallocate (ky_phiba)
! *******************
!     
!     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)
      
    wtif(ne,:,:)=0.5*(den2m(:,:)+v2m(:,:)+0.5*tpar2m(:,:)+tp2m(:,:))
    wkif(ne,:,:)=v2m(:,:)
    wpif(ne,:,:)=pr2m(:,:)
!     
!  Phi* Phi    for each m and n at z=0
!     
    l=ldb/2
    psp(ne,:,:)=real(potential(l,:,:))**2+aimag(potential(l,:,:))**2
!     
! Other variables we want to save     
!     
    timo(ne)=time(md)
!
! Flux surface averaged poloidal and toroidal flows, for m=0,n=0 mode
!
    if(iperiod /= 2) then
       l_left=1
       l_right=ld-1
    endif
  
    if(.not.associated(dz)) allocate (dz(l_left:l_right))
    dz0=2.*x0/ldb 
    do l=l_left,l_right
!     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
! rkx, bz not upgraded for general geometry yet, so skip this part
!     
    if(igeo == 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)*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 
             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/qsf
             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
!    write(*,*) 've/ve0: ',veh0(ne)/veh0(1)
!
!     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
!     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)
             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
             mr=mrr(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
       if(.not.associated(ninv)) allocate (ninv(-nd:nd))
       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(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 <= 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 >= -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 
  
    allocate (phit0(md, nd))
    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
    deallocate (phit0)
    
    init=0
    
    deallocate (dummy, dum1k, dum2k, den2m, &
         v2m, tpar2m, tp2m, pr2m, d1eta2, d2eta2, d3eta2, d1par2, d2par2, &
         wupsi2, d1r2, d3r2, d4r2, d5r2, d6r2, d7r2, d8r2, d9r2)

    deallocate (totpr, ky_psiv, nbar, &
         ky_psiq, p_perp, kyn1, kyt1, kyt2, kwn1, kwt1, kwt2, &
         kpt1, kpt2, ktu1, ktt1, ktt2, ktr1, ktt3, ktt4)

  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
!
 
    real, allocatable, dimension(:,:) :: dummy, u2
    real dum
!    real nuur,nuui
    integer nea
    
    real den2,v2,tpar2,tp2
    real d1eta,d2eta,d3eta,d1par,d2par
    real d1r,d3r,d4r,d5r,d6r,d7r,d8r,d9r
!
    allocate (dummy(md, nd), u2(md,nd))
    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 
    grtmx(nea,:,:) = u2(:,:)
    
    pwenr=0.5*(den2+v2+tpar2/2.+tp2)
    
    deallocate (dummy, u2)

  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.
!     
!
    real, allocatable, dimension(:,:) :: u2, dummy
    
    real edif,etot,dum
    real den2,v2,tpar2,tp2,wenrn,dedt,dt1
    integer m,n 
    
    allocate (u2(md, nd))
    call volsq(potential,potential,dum,u2)
    
    if(ncycle > 2) then 
       do n=1,nd 
          do m=1,md
             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 
       grtmx(ne,:,:)=0.0 
    endif
    deallocate (u2)
!
!     calculate energy conservation, 2nd order accurate in dt
!     
    allocate (dummy(md, nd))
    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)
    deallocate (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:
!
    dt1=dt(md)
    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!
!     

    complex, dimension(:,:,:) :: a, b
    real, dimension(:,:) :: axb_by_mode
    real, dimension(:), pointer, save :: dz, zero
    
!    integer :: init = 1
    real axb,vnorm,dz0
    integer l,m,n
    
    if(.not.associated(dz)) allocate (dz(l_left:l_right), zero(md))
    dz0=2.*x0/ldb
    if(iperiod /= 2) then
       l_left=1
       l_right=ld-1
    endif
    
    axb=0.0
!      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
       if(igeo == 0) then
          dz(l)=dz0
       else
          dz(l)=dz0*jacobian(l)
       endif
    enddo
    
    do n=1,nd
       do m=1,md
          axb_by_mode(m,n)=0.
          vnorm=0.
          do l=l_left,l_right
             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
!
! 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.)

    complex, dimension(:,:,:) :: a, b
    real, dimension(:,:) :: axb_by_mode
    real, dimension(:), pointer, save :: dz, zero
    
    real, dimension(:,:,:), allocatable :: c, d
    
! old coding had a bug that probably never showed up since kz = lz
!      real axb,c(kz,mz,nz),d(lz,mz,nz)
    real axb, vnorm
    integer l,m,n,j,l0
    
    if(.not.associated(dz)) allocate (dz(l_left:l_right), zero(md))
    if(iperiod /= 2) then
       l_left=1
       l_right=ld-1
    endif
    
    do l=l_left,l_right
       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
!
    allocate (c(kdpass, md, nd), d(kdpass, md, nd))
    do n=1,nd
       do m=2,md
          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
       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
    enddo
!
!     now average in kappa
!
    do m=2,md
       do n=1,nd
          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=1,nd
       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
    enddo
!
!     finally volume average
!
    axb=0.0
    do n=1,nd
       do m=1,md
          axb_by_mode(m,n)=0.
          vnorm=0.
          do l=l_left,l_right
             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
    
    deallocate (c, d)
    
  end subroutine volsqk
  
  subroutine volflux(a, b, axb, axb_by_mode, nsp)
!
! 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.  
!     

    complex, dimension(:,:,:,:) :: a, b
    real, dimension(:,:) :: axb_by_mode
    real, dimension(:), pointer, save :: dz, zero
    real :: axb, vnorm
    integer l, m, n, nsp

    if(.not.associated(dz)) allocate (dz(l_left:l_right), zero(md))
!
! This logic is wrong for input md=-1.     
!
    zero(1)=1.0
    do m=2,md
       zero(m)=0.5
    enddo
    
    if(iperiod /= 2) then
       l_left=1
       l_right=ld-1
    endif
  
    do l=l_left,l_right
       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
    do n=1,nd
       do m=1,md
          axb_by_mode(m,n)=0.
          vnorm=0.
          do l=l_left,l_right
             vnorm=vnorm+dz(l)
             axb_by_mode(m,n)=axb_by_mode(m,n) &
                  +(real(a(l,m,n,nsp))*real(b(l,m,n,nsp)) &
                  +aimag(a(l,m,n,nsp))*aimag(b(l,m,n,nsp)))*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

  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.)

    complex, dimension(:,:,:) :: a, b
    real, dimension(:,:) :: axb_by_mode
    real, dimension(:), pointer, save :: dz, zero
    real, allocatable, dimension(:,:,:) :: c, d
    
    real vnorm, axb
    integer l,m,n,j,l0
    
    if(.not.associated(dz)) allocate (dz(l_left:l_right), zero(md))
    if(iperiod /= 2) then
       l_left=1
       l_right=ld-1
    endif
    
    do l=l_left,l_right
       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)=2.0
       zero(m)=0.5
    enddo
!
!     first multiply pitch angle functions together
!
    allocate (c(kdpass, md, nd), d(kdpass, md, nd))
    do n=1,nd
       do m=2,md
          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
       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
    enddo
!
!     now average in kappa
!
    do m=2,md
       do n=1,nd
          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=1,nd
       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
    enddo
!
!     finally volume average
!
    axb=0.0
    do n=1,nd
       do m=1,md
          axb_by_mode(m,n)=0.
          vnorm=0.
          do l=l_left,l_right
             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

    deallocate (c, d)

  end subroutine volkflux

end module diagnostics