c INCLUDE file itg.cmn:  ITG subroutine common blocks
c
c     (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, 
c     and Gregory W. Hammett. ALL RIGHTS RESERVED.
c
c major parameters for the grid, fields, and modes
c
      integer ld	! # of parallel meshpoints
      integer ldb	! (# of parallel mespoints) - 1
      integer kd	! # of kappa grid points (pitch angle) trapped	
      integer kdpass	! # of kappa grid points (pitch angle) +passing
      integer nffty	!  # of y grid points, with dealiasing
      integer nfftx	!  # of x grid points, with dealiasing
      integer nmin	! fourier modes range from nmin to nmax.
      integer nmax      ! nmin=nmax=0 for single helicity
      integer nd	! total # of n's,  nd=nmax-nmin+1
      integer nspecies  ! number of species
      integer iperiod   ! 1 -> periodic in theta b.c.'s, 0 -> 0 b.c.'s
			! 2 -> extensions, 3 -> connected FFTs
      integer igeo	! switch for general geometry
      integer iflow     ! switch for treatment of m=0 modes
      integer icrit     ! switch for L_Tcrit, D_ML searches
      integer idt	! switch for adjustable time step
      integer jrad      ! flux surface label (for multiple calls to itgc)
      real pi		! 3.14159...
      real dtadj	! dt -> 1.05*dtadj*dt or dt -> dt/dtadj
      real dt0		! time step (in units of Ln/cs)
      integer lin	! 1= linear run, 0=nonlinear run.

      integer n_p          ! number of particles per cloud
      integer n_zpic       ! number of theta grid points in PIC code
      real bmax            ! bmax of PIC code
      real chi_int 	   ! int(chi dt)
      logical binary       ! true -> binary i/o for big files

      character layout*6 ! determines layout of k_r 
      character runname*80 ! Name of this run.  Associated i/o files
                        ! will be named RUNNAME.in, etc.
      character fullname*80 ! full path name of input file
      integer lrunname    ! length of runname
      character note*80 ! note about the run.

      integer lhistory	! phi @ r(lhistory) saved for plots
	 
      real etae		! eta_e
      real rmime	! m_i / m_e
      real r(lz)	! grid for the radial (x) coordinate
      real kap(kz+1)	! grid for kappa (pitch angle) coordinate
      real dkap(kz)     ! dkappa for bounce averaging
      real taub(kz)	! bounce time, trapped and passing
      real pmat(lz,lz)  ! matrix for poisson inversion	
c 
c For electron collisions:
c
      real bratio(kz)	! bounce ave. of (B-Bmin)/B
      real jac(kz)	! Jacobian
      real diff(0:kz)
      real aec(kz,mz),bec(kz,mz),cec(kz,mz) ! for electron collision operator 
      real zbeta(kz,mz),zgamma(kz,mz)	    ! for electron collision operator 
c
      real bdel(lz)	! By/B profile (includes shear effects)
      real omegad(lz,mz,nz,nspecz) 	! omega_d profile (toroidal terms)
      real omegade(kz,mz,nz) 	! omega_de profile (electron terms)
      real b_mag(lz)	! magnitude of B profile
      real bmaginv(lz)	! 1/B profile
      real bgrad(lz)    ! 1/B dB/dl profile
      real kpar(lz*nconz,ndomz)	! array for kpar fft's
      real stochf	! stochasticity parameter for poisson	
c
c mdspec doesn't need to be in a common block? gwh
	integer mmin(nz),mmax(nz),md,mdspec,nddp,nddm,ivex(lz,mzz),
     *  malias,nalias,madd,nadd,mlow,mhi,iseed,n0
c
c FOURIER TRANSFORM CONVENTIONS FOR THE ITG CODE
c
c ?? The description below is for the old slab version of the code,
c and needs updating for the toroidal version:
c
C The code solves fluid equations using a spectral representation for
C the poloidal (y) and toroidal (z) directions, and a second-order
C finite difference method for the radial (x) direction.  All
C quantities can be expanded in Fourier harmonics as:
c
c	wfield(x,y,z) = Sum_mn w(x,m,n) exp( + i m y / y0 - i n z / z0)
c
c where w(x,m,n) are complex coefficients.  
c
c The (-m,-n) mode is linearly dependent on the (m,n) mode, so
c we can eliminate half of the modes.  We choose to keep just
c the modes with m >= 0, although n <0, =0, and >0 are all be kept.
c Caution: though the m>0 modes are all independent, the m=0 modes 
c are not.
c
c All primary variables are for guiding centers, normalized to the
c equilibrium value for each species, i.e. n(j)=(L_ne/rho_i) n_j1/n_j0, 
c except the potential, which is Phi=(L_ne/rho_i) e Phi_1/T_i0, where 
c T_i0 is the main species equilibrium temperature, rho_i=v_ti/Omega_i,
c and v_ti=sqrt(T_i0/m_i)  (E.g., 'density' is the guiding center density, 
c which is different from the particle density due to FLR effects.)  
c
c Time is normalized to L_ne/v_ti. 
c
c 	density  = guiding center density n
c 	u_par  = parallel velocity
c 	t_par = parallel temperature
c 	q_par = parallel flux of parallel heat
c 	t_perp  = perpendicular temperature
c 	q_perp  = parallel flux of perpendicular heat
c 	potential  = phi (electrostatic potential)
c
c The "home" basis for most fields is (theta,ky,kx). For example,
c density(lz,mz,nz,nspecz) is the guiding center density vs.
c     l=1,ldb  along-the-field-line grid index
c     m=1,md   ky (poloidal mode) index
c     n=1,nd   kx (radial mode) index
c     nsp=1,nspecies  species index

      complex density(lz,mz,nz,nspecz),u_par(lz,mz,nz,nspecz)
      complex t_par(lz,mz,nz,nspecz),q_par(lz,mz,nz,nspecz)
      complex t_perp(lz,mz,nz,nspecz),q_perp(lz,mz,nz,nspecz)
      complex potential(lz,mz,nz)
      complex phiav_old(mz,nz) 	
c
c Electron variables in pitch angle space
c
      complex e_density(kz,mz,nz),e_p(kz,mz,nz),e_r(kz,mz,nz)
      complex e_t(kz,mz,nz)
      complex phi_ba(kz,mz,nz)
c
c Electron variables in theta space
c
      complex phi_bk(lz,mz,nz),e_denk(lz,mz,nz)
c
c Nonlinear terms.  Could be condensed down to one per equation for 
c a small memory cost.
c
      complex nl_density1(lz,mz,nz,nspecz),nl_density2(lz,mz,nz,nspecz),
     . nl_upar1(lz,mz,nz,nspecz),nl_upar2(lz,mz,nz,nspecz),
     . nl_tpar(lz,mz,nz,nspecz),nl_qpar(lz,mz,nz,nspecz),
     . nl_tperp1(lz,mz,nz,nspecz),nl_tperp2(lz,mz,nz,nspecz),
     . nl_tperp3(lz,mz,nz,nspecz),nl_qperp1(lz,mz,nz,nspecz),
     . nl_qperp2(lz,mz,nz,nspecz),nl_qperp3(lz,mz,nz,nspecz)
c
c Electron nonlinear terms in pitch angle space
c
      complex nl_edensity(kz,mz,nz),nl_ep(kz,mz,nz),nl_er(kz,mz,nz)
      complex nl_et(kz,mz,nz)
c
c The following quantities include FLR corrections to Phi:
c phi_n is the FLR-corrected w_* term for the density equation
c phi_u is the FLR-corrected w_* term for the parallel velocity equation
c phi_tperp is the FLR-corrected w_* term for the perp. temperature equation
c phi_qperp is the FLR-corrected w_* term  for the par-perp heat equation
c phi_nd is the FLR-corrected w_d term for the density equation
c phi_tpard is the FLR-corrected w_d term for the perp. temperature equation
c phi_tperpd is the FLR-corrected w_d term  for the par temp. equation
c etc.
c
      complex phi_n(lz,mz,nz,nspecz),phi_u(lz,mz,nz,nspecz),
     . phi_tpar(lz,mz,nz,nspecz),phi_tperp(lz,mz,nz,nspecz),
     . phi_qperp(lz,mz,nz,nspecz),phi_rperp(lz,mz,nz,nspecz),
     . phi_flr1(lz,mz,nz,nspecz),phi_flr2(lz,mz,nz,nspecz),
     . phi_flr3(lz,mz,nz,nspecz),phi_nd(lz,mz,nz,nspecz),
     . phi_tpard(lz,mz,nz,nspecz),phi_tperpd(lz,mz,nz,nspecz),
     . phi_qperpb(lz,mz,nz,nspecz)

      real lambda1,lambda2,nu1,nu2,cflx,cfly,vmax,rdiff_0
      real vy,vz,shr,qsf,epsn,eps,epse,alpha,gamma_0,rdiff(mz,nz)
      real etai,nuii,nueeff,rmu1,rootm,rootn,nu_ii(nspecz),tol
      real tiovte,rkkin,rkkperp,etaipar,rkkperp2,
     * nuar,nuai,nubr,nubi,nucr,nuci,nudr,nudi,
     * nu1r,nu1i,nu2r,nu2i,nu3r,nu3i,nu4r,nu4i,nu5r,nu5i,nu6r,nu6i, 
     * nu7r,nu7i,nu8r,nu8i,nu9r,nu9i,nu10r,nu10i,nu11r,nu11i,
     * nu12r,nu12i,nu13r,nu13i,nu14r,nu14i,zDq1,zBq1,zEq1,s_par,
     * s_perp
      integer mrr(mz,nz),nrr(nz),igradon,iphi00,iflr,nstp,ifilter,ihdf,
     * nprnt,nread,movieon,ninterv,meq,neq,
     * navgflag,ncycle,ne,nfreqcnt,nperpmom,nparmom,nemom,iodd,ntrace,
     * iexp,inlpm,l_left,l_right,stepnlps
      integer ikx ! Obsolete.  Used to control pseudo-spectral vs.
                  ! finite-difference d/dx in nonlinear terms.
      real dr(lz),r0(mz,nz),x0,xp,y0,z0,pmag,dens,avgflag,vxmax,vymax

      real timo(nez),wpfx(nez,nspecz),gamx(nez),wrhx(nez),wrhy(nez)
      real fluxe(nez),fluxi(nez,nspecz),qfluxe(nez)
      real wenr,wenr1,wiox(nez,mz,nz),phik2(lz,nez),pwenr,dtold,drive,
     * pcerr(nez),mgamx(mz,nz,nez)
      real therm(lz),wxsp(nez),grtmx(nez,mz,nz),
     * wplx(nez,mz,nz),wakx(nez),waky(nez),
     * wkif(nez,mz,nz),wpif(nez,mz,nz),wtif(nez,mz,nz),
     * wenx(nez),wpcx(nez),wppif(nez,mz,nz),
     * wlux(nez),wdpx(nez),wdvx(nez),wdfx(nez),
     * wdux(nez),wu2x(nez),wp2x(nez),wv2x(nez),
     * timf(nez),upol(nez),utor(nez),
     * psp(nez,mz,nz),phirms(nez)
      real phisq(nez,mz,nz)		! potential^2 for each mode
      real fluximn(nez,mz,nz,nspecz) 	! particle flux for each mode
      real qfluximn(nez,mz,nz,nspecz)	! heat flux for each mode
      real fluxemn(nez,mz,nz)		! electron particle flux for each mode
      real qfluxemn(nez,mz,nz)		! electron heat flux for each mode
      real cx(nez,4*nz),cy(nez,4*mz),cz(nez,lz)
      real czc(nez,lz),czcn(nez,lz)
      real czn(nez,lz),cznn(nez,lz)
      complex phi00(nez,nz),den00(nez,nz),upar00(nez,nz)
      complex tpar00(nez,nz),tperp00(nez,nz)
      real ct(nez)
      real drt(14,nez)	
      real wenrk(nez,mz,nz),dketa(nez,mz,nz),dkpar(nez,mz,nz)
      real dktor(nez,mz,nz),dktdp(nez,mz,nz),wkups(nez,mz,nz)
      real dktot(nez,mz,nz)

      complex utim(nez,mz,nz),phit0(mz,nz)
      integer mfr,nfr,nfreq,ntotal,ntimf
      character*30 date

      real rky(mz,nz)       		! ky for each mode
      real rky2(mz,nz)    		! ky**2 for each mode.
      real rkx2(lz,mz,nz)		! kx**2 grid.
      real rkx(lz,mz,nz)		! kx grid.
      real rkperp2(lz,mz,nz)		! k_perp**2 

      real gnwgt(lz,mz,nz,nspecz)	! density equation weight
      real gtpwgt(lz,mz,nz,nspecz)	! T_par equation weight
      real gtwgt(lz,mz,nz,nspecz)	! T_perp equation weight
      real gnwgtd(lz,mz,nz,nspecz)	! density equation weight, toroidal 
      real gtpwgtd(lz,mz,nz,nspecz)	! T_par equation weight, toroidal 
      real gtwgtd(lz,mz,nz,nspecz)	! T_perp equation weight, toroidal 
      real gqwgtb(lz,mz,nz,nspecz)	! q_perp equation weight, toroidal 
      real grwgt(lz,mz,nz,nspecz)	! r_perp equation weight
      real g0wgt(lz,mz,nz,nspecz)	! gamma_0^(1/2)
      real g0wgtp(lz,mz,nz,nspecz)	! gamma_0^(1/2) in poisson
      real flrwgt(lz,mz,nz,nspecz)	! nonlinear FLR weight
      real flrwgt2(lz,mz,nz,nspecz)	! T nonlinear FLR weight
      real flrwgt3(lz,mz,nz,nspecz)	! n,u nonlinear FLR weight

      real gewgt(lz,mz,nz,nspecz)
      real twgt(lz,mz,nz,nspecz)
      real nwgt(lz,mz,nz,nspecz)
      real pfilter(lz,mz,nz)		! filter to get phi from n.
      real pfilter2(lz,mz,nz)		! filter to get phi from n.
      real pfilter3(lz,mz,nz,nspecz)	! filter to get pol. density from n
      real denom0(mz,nz)		! field-line-average term

      real wgtba1(kz,lz)		! weight for bounce avg.	
      real wgtba2(kz,lz)		! weight for kappa avg.	

c
c Geometrical coefficients
c
      real gbd(lz),gbd0(lz),cvd(lz),cvd0(lz),gpar,jacobian(lz)
      real gd2(lz),gd21(lz),gd22(lz),gradrho(lz)

      character*10 initx,initkx,inity,initky,initkp
c      real trigx(2*nzz) 	    	! NAG workspace
c      real trigkx(2*nzz)	    	! NAG workspace
c      real trigy(2*mzz) 	    	! NAG workspace
c      real trigky(2*mzz)	    	! NAG workspace
c      real trigkp(2*lz)		    	! NAG workspace
      real tablekp(ndomz,101+2*lz*nconz)       	! mcfft trig array
      real tablex(100+2*nzz)        	! mcfft trig array		
      real tabley(100+2*mzz)        	! mcfft trig array		
      integer ncon(ndomz),scon(ndomz),conpos(mz,nz),mcon(mz,nz)
      integer icon(mz,nz),ndom,mmcon(ndomz,mz*nz),nncon(ndomz,mz*nz)
	

c  Added by GWH
      real worky((2*mzz+4)*lz*nzz)  ! fft workspace

c	ion nonlinear storage arrays
      real kxa3(mzz,nzz*lz),kya3(mzz,nzz*lz)
      real kxb3(mzz,nzz*lz),kyb3(mzz,nzz*lz)
c	electon nonlinear storage arrays
      real kxae3(mzz,nzz*kz),kyae3(mzz,nzz*kz)
      real kxbe3(mzz,nzz*kz),kybe3(mzz,nzz*kz)

      real rho(nspecz),tau(nspecz),charge(nspecz),n_I(nspecz),
     . rmass(nspecz),eta(nspecz),eta_par(nspecz),vt(nspecz),
     . Ln(nspecz)

      real etak(mz,nz,nspecz),etak_par(mz,nz,nspecz),dm_out(mz,nz)

      common/geom/gbd,gbd0,cvd,cvd0,gpar,gd2,gd21,gd22,jacobian,
     .	gradrho
      common/dmix1/etak,etak_par,gamma_0,rdiff_0,rdiff,tol
      common/dmix2/dm_out
      common/conn/ncon,scon,conpos,mcon,icon,ndom,mmcon,nncon	
      common/multispecies/rho,tau,charge,n_I,rmass,eta,eta_par,
     .	vt,Ln,nu_ii
      common/nlpssave/kxa3,kya3,kxb3,kyb3,vxmax,vymax,stepnlps
      common/nlpsesave/kxae3,kyae3,kxbe3,kybe3	
c      common/ffttrig/trigx,trigkx,trigy,trigky,trigkp
      common/mcfft/tablekp,tablex,tabley,worky
      common/fftinit/initx,initkx,inity,initky,initkp
      common/flrcom/rky,rky2,rkx2,rkperp2,rkx,gnwgt,gtpwgt,gtwgt,grwgt,
     * g0wgt,flrwgt,flrwgt2,flrwgt3,gewgt,twgt,nwgt,pfilter,pfilter2,
     * g0wgtp,denom0,gnwgtd,gtwgtd,gtpwgtd,gqwgtb,pfilter3
      common/moments/density,u_par,t_par,q_par,t_perp,q_perp,potential
      common/electrons/e_density,e_p,e_r,e_t,phi_ba,e_denk,phi_bk	
      common/nonlinear/nl_density1,nl_density2,nl_upar1,nl_upar2,
     * nl_tpar,nl_qpar,nl_tperp1,nl_tperp2,nl_tperp3,nl_qperp1,
     * nl_qperp2,nl_qperp3
      common/nonlineare/nl_edensity,nl_ep,nl_er,nl_et
      common/flr_terms/phi_n,phi_u,phi_tpar,phi_tperp,phi_qperp,
     * phi_rperp,phi_flr1,phi_flr2,phi_flr3,phi_nd,phi_tpard,phi_tperpd,
     * phiav_old,phi_qperpb
      common /ddat/ld,ldb,kd,kdpass,nffty,nfftx,nmin,nmax,nd,lin,
     . lhistory,
     . 	iexp,nddp,nddm,iseed,n0,jrad,n_p,n_zpic,bmax,chi_int,binary
      common /ddat3/pi,dt0,tiovte,pcerr,etaipar,rootn,rootm,dtadj
      common /ddat2/mmin,mmax,md,ifilter,ikx,mdspec,idt,inlpm,icrit,
     .	igeo,iflow
      common /dpara/vy,vz,shr,qsf,epsn,eps,epse,alpha,
     * etai,nuii,nueeff,rmu1,mlow,mhi,
     * rkkin,rkkperp,malias,nalias,nadd,madd,rkkperp2,
     * lambda1,lambda2,nu1,nu2,cflx,cfly,vmax,
     * nuar,nuai,nubr,nubi,nucr,nuci,nudr,nudi,
     * nu1r,nu1i,nu2r,nu2i,nu3r,nu3i,nu4r,nu4i,nu5r,nu5i,nu6r,nu6i, 
     * nu7r,nu7i,nu8r,nu8i,nu9r,nu9i,nu10r,nu10i,nu11r,nu11i,
     * nu12r,nu12i,nu13r,nu13i,nu14r,nu14i
      common /grid/r,kap,dkap,bdel,omegad,omegade,b_mag,bmaginv,bgrad,
     * kpar,mrr,nrr,l_left,l_right,dr,r0,x0,xp,y0,z0,stochf,
     * rmime,etae
      common /bounce/ taub,wgtba1,wgtba2,pmat
      common/ecoll/bratio,jac,diff,aec,bec,cec,zbeta,zgamma
      common /mesh/zDq1,zBq1,zEq1,s_par,s_perp
      common /control/avgflag,pmag,dens,igradon,iphi00,iflr,ihdf,
     * nprnt,nread,movieon,ntrace,nstp,nperpmom,nemom,
     * ninterv,meq,neq,navgflag,ncycle,ne,nfreqcnt,nparmom,iodd,
     * nspecies,ivex,iperiod,ntimf
      common /pdata0/wenr,wenr1,wiox,phik2,timo,wpfx,gamx,wrhx,
     * fluxe,fluxi,qfluxe,fluximn,qfluximn,phisq,fluxemn,qfluxemn
      common /pdata1/wrhy,mgamx,wdux,wu2x,wp2x,wv2x,wlux,wdpx,
     * wdvx,wdfx,
     * wenx,wpcx,wppif,wkif,wpif,wtif,wplx,wakx,waky,therm,wxsp
      common /pdata2/ grtmx,psp,phirms,cx,cy,cz,czc,czcn,czn,cznn,ct,
     * phi00,den00,upar00,tpar00,tperp00,drt,pwenr,dtold,drive
      common /pdata3/ wenrk,dketa,dkpar,dktor,dktdp,wkups,dktot,phit0,
     * utim,mfr,nfr,timf,nfreq,ntotal
      common /flow/utor,upol

c When doing debug runs, it is sometimes useful to force the plots
c to have the same date/time label as a previous run, so that one
c can verify results by typing "diff run1.m run2.m".  In order to
c do this, set debug_plotlabel to the previous date/time label:
      character debug_plotlabel*30

c character and non-character types can't be mixed in common blocks
c on the CM:
      common /strings/ note, runname, fullname, date, debug_plotlabel
      common /strings/ layout
      common /stringsl/ lrunname