module itg_data
!
!     (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, 
!     and Gregory W. Hammett. ALL RIGHTS RESERVED.


  complex, parameter :: ii = (0., 1.)

!
! major parameters for the grid, fields, and modes
!
  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 ndomz      ! # of different FFT lengths.  
  integer nconz      ! # of connected FFT domains
  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 :: md, nddp, nddm, malias, nalias, madd, nadd, mlow, mhi, iseed, n0

  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
  real stochf	     ! stochasticity parameter for poisson

  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
  integer :: mfr, nfr, nfreq, ntotal, ntimf
  character*30 date
	 
  real :: etae		           ! eta_e
  real :: rmime	                   ! m_i / m_e

  real, dimension(:), allocatable :: time, dt
  
  real, dimension(:), allocatable :: r, dr! grid for the radial (x) coordinate 
  real, dimension(:), allocatable :: kap  ! grid for kappa (pitch angle) coordinate 
  real, dimension(:), allocatable :: dkap ! dkappa for bounce averaging
  real, dimension(:), allocatable :: taub ! bounce time, trapped and passing
!  real, dimension(:,:), allocatable :: pmat ! matrix for poisson inversion
! 
! For electron collisions:
!
! bratio:       bounce avg. of (B-Bmin)/B
! jac:          electron Jacobian
! diff:         parameter for bounce averaging operator

  real, dimension(:), allocatable :: bratio, jac, diff

! omegad:       omega_d vs. theta (toroidal terms)
! omegade:      omega_de vs. theta (electron terms)
! b_mag:        |B| vs. theta
! bmaginv:      1/B vs. theta
! bgrad:        1/B dB/dl vs. theta
! kpar:         array for k_parallel FFT's
  real, dimension(:,:,:,:), allocatable :: omegad 
  real, dimension(:,:,:),   allocatable :: omegade     
  real, dimension(:,:),     allocatable :: kpar     
  real, dimension(:),       allocatable :: b_mag, bmaginv, bgrad

! Many versions ago, a wedge of modes was kept; mmin, mmax specified the wedge.
! nrr and mrr contain the actual n and m numbers for a given index value.
!
  integer, dimension(:), allocatable ::  mmin, mmax
  integer, dimension(:), allocatable :: nrr
  integer, dimension(:,:), allocatable :: mrr

!
! FOURIER TRANSFORM CONVENTIONS FOR THE ITG CODE
!
! ?? The description below is for the old slab version of the code,
! and needs updating for the toroidal version:
!
! The code solves fluid equations using a spectral representation for
! the poloidal (y) and toroidal (z) directions, and a second-order
! finite difference method for the radial (x) direction.  All
! quantities can be expanded in Fourier harmonics as:
!
!	wfield(x,y,z) = Sum_mn w(x,m,n) exp( + i m y / y0 - i n z / z0)
!
! where w(x,m,n) are complex coefficients.  
!
! The (-m,-n) mode is linearly dependent on the (m,n) mode, so
! we can eliminate half of the modes.  We choose to keep just
! the modes with m >= 0, although n <0, =0, and >0 are all be kept.
! Caution: though the m>0 modes are all independent, the m=0 modes 
! are not.
!
! All primary variables are for guiding centers, normalized to the
! equilibrium value for each species, i.e. n(j)=(L_ne/rho_i) n_j1/n_j0, 
! except the potential, which is Phi=(L_ne/rho_i) e Phi_1/T_i0, where 
! T_i0 is the main species equilibrium temperature, rho_i=v_ti/Omega_i,
! and v_ti=sqrt(T_i0/m_i)  (E.g., 'density' is the guiding center density, 
! which is different from the particle density due to FLR effects.)  
!
! Time is normalized to L_ne/v_ti. 
!
! 	density  = guiding center density n
! 	u_par  = parallel velocity
! 	t_par = parallel temperature
! 	q_par = parallel flux of parallel heat
! 	t_perp  = perpendicular temperature
! 	q_perp  = parallel flux of perpendicular heat
! 	potential  = phi (electrostatic potential)
!
! The "home" basis for most fields is (theta,ky,kx). For example,
! density(lz,mz,nz,nspecz) is the guiding center density vs.
!     l=1,ldb  along-the-field-line grid index
!     m=1,md   ky (poloidal mode) index
!     n=1,nd   kx (radial mode) index
!     nsp=1,nspecies  species index

  complex, dimension(:,:,:,:), allocatable :: density, u_par, t_par, q_par, t_perp, q_perp
  complex, dimension(:,:,:),   allocatable :: potential
  complex, dimension(:,:),     allocatable :: phiav_old

!
! Electron variables in pitch angle space
!
  complex, dimension(:,:,:), allocatable :: e_density, e_p, e_r, e_t, phi_ba

!
! Electron variables in theta space
!
  complex, dimension(:,:,:), allocatable :: phi_bk, e_denk
!
! Nonlinear terms.  Could be condensed down to one per equation for 
! a small memory cost.
!
  complex, dimension(:,:,:,:), allocatable :: nl_density1, nl_density2, &
       nl_upar1, nl_upar2, nl_tpar, nl_qpar, nl_tperp1, nl_tperp2, &
       nl_tperp3, nl_qperp1, nl_qperp2, nl_qperp3
!
! Electron nonlinear terms in pitch angle space
!
  complex, dimension(:,:,:), allocatable :: nl_edensity, nl_ep, nl_er, nl_et

  real :: lambda1, lambda2, nu1, nu2, cflx, cfly, vmax, rdiff_0
  real :: vy, vz, shr, qsf, epsn, eps, epse, alpha, gamma_0
  real, dimension(:,:), allocatable :: rdiff

  real :: etai, nuii, nueeff, rmu1, rootm, rootn, 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 :: igradon, iphi00, iflr, nstp, ifilter, ihdf, nprnt, nread, &
       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, dimension(:,:), allocatable :: r0
  real :: x0, xp, y0, z0, pmag, dens, avgflag, vxmax, vymax

  real, dimension(:), allocatable :: timo, gamx, wrhx, wrhy, &
       fluxe, qfluxe, pcerr, wxsp, wakx, waky, wenx, wlux, &
       wdpx, wdvx, timf, upol, utor, phirms

  complex, dimension(:), allocatable :: phih0, uparh0, denh0, &
       tparh0, tperph0, qparh0, qperph0, veh0, uparph0

  real, dimension(:,:),   allocatable :: wpfx, fluxi
  real, dimension(:,:,:), allocatable :: grtmx, wkif, wpif, wtif, psp, phisq
  real, dimension(:,:,:), allocatable :: mgamx

  real :: wenr, wenr1, pwenr, dtold, drive

  ! ion particle and heat flux for each mode: 
  real, dimension(:,:,:,:), allocatable :: fluximn, qfluximn 
  ! electron particle and heat flux for each mode: 
  real, dimension(:,:,:), allocatable :: fluxemn, qfluxemn

! correlation function variables
  real, dimension(:,:), allocatable :: cx, cy, cz, czc, czcn, czn, cznn

  complex, dimension(:,:), allocatable :: phi00, den00, upar00, tpar00, tperp00

  real, dimension(:),     allocatable :: ct
  real, dimension(:,:),   allocatable :: drt
  real, dimension(:,:,:), allocatable :: wenrk, dketa, dkpar, dktor, dktdp, wkups, dktot

  complex, dimension(:,:,:), allocatable :: utim

  real, dimension(:,:),   allocatable :: rky, rky2  ! ky, ky**2 for each mode
  real, dimension(:,:,:), allocatable :: rkx, rkx2, rkperp2 ! kr, kr**2, kperp**2 for each mode

  ! Weights related to FLR effects
  
  real, dimension(:,:,:,:), allocatable :: gnwgt, gtpwgt, gtwgt, gnwgtd, gtpwgtd, gtwgtd, &       
       gqwgtb, grwgt, g0wgt, g0wgtp, flrwgt, flrwgt2, flrwgt3, gewgt, twgt, nwgt, pfilter3

  real, dimension(:,:,:), allocatable :: pfilter

  real, dimension(:,:), allocatable :: denom0	! field-line-average term
  real, dimension(:,:), allocatable :: wgtba1	! weight for bounce avg.
  real, dimension(:,:), allocatable :: wgtba2	! weight for kappa avg.	

!
! Geometrical coefficients
!
  real :: gpar
  real, dimension(:), allocatable :: gbd, gbd0, cvd, cvd0, jacobian, gd2, gd21, gd22, gradrho     

  character*10 initx,initkx,inity,initky,initkp

  real, dimension(:,:), allocatable :: tablekp    ! mcfft trig array

  integer, dimension(:,:), allocatable :: conpos, mcon, icon, nncon, mmcon
  integer, dimension(:),   allocatable :: ncon, scon
  integer ndom
	
  real, dimension(10) :: rho, tau, charge, n_I, rmass, eta, eta_par, vt, Ln, nu_ii       
  real, dimension(:,:,:), allocatable :: etak, etak_par
  real, dimension(:,:),   allocatable :: dm_out

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

  integer :: mplot2, nplot2    ! (m,n) mode number to make plots for
  integer, dimension(:), allocatable :: dash ! dash pattern to identify mode # in plots
  
contains

  subroutine allocate
!
! Allocate most arrays for gryffin.  
!
    allocate (density(ld, md, nd, nspecies),  u_par(ld, md, nd, nspecies),  &
                t_par(ld, md, nd, nspecies),  q_par(ld, md, nd, nspecies),  &
               t_perp(ld, md, nd, nspecies), q_perp(ld, md, nd, nspecies),  &
               omegad(ld, md, nd, nspecies))
    
    allocate (potential(ld, md, nd),  phiav_old(md, nd))

    allocate (nl_density1(ld, md, nd, nspecies),  nl_density2(ld, md, nd, nspecies),&
          nl_upar1(ld, md, nd, nspecies), nl_upar2(ld, md, nd, nspecies),  &
           nl_tpar(ld, md, nd, nspecies),  nl_qpar(ld, md, nd, nspecies),  &
         nl_tperp1(ld, md, nd, nspecies), nl_tperp2(ld, md, nd, nspecies), &
         nl_tperp3(ld, md, nd, nspecies), nl_qperp1(ld, md, nd, nspecies), &
         nl_qperp2(ld, md, nd, nspecies), nl_qperp3(ld, md, nd, nspecies))

    if(epse > 0.) then
       allocate (e_density(kdpass, md, nd),  e_p(kdpass, md, nd),  &
                       e_r(kdpass, md, nd),  e_t(kdpass, md, nd),  &
                    phi_ba(kdpass, md, nd))

       allocate (nl_edensity(kd, md, nd), nl_ep(kd, md, nd), &
                       nl_er(kd, md, nd), nl_et(kd, md, nd))
    endif

    allocate (phi_bk(ld, md, nd),  e_denk(ld, md, nd))
!    allocate (pmat(ld,  ld))  ! not used

    allocate (gnwgt(ld, md, nd, nspecies),  gtpwgt(ld, md, nd, nspecies),  &
              gtwgt(ld, md, nd, nspecies),  gnwgtd(ld, md, nd, nspecies),  &
            gtpwgtd(ld, md, nd, nspecies),  gtwgtd(ld, md, nd, nspecies),  &       
             gqwgtb(ld, md, nd, nspecies),   grwgt(ld, md, nd, nspecies),  &
              g0wgt(ld, md, nd, nspecies),  g0wgtp(ld, md, nd, nspecies),  &
             flrwgt(ld, md, nd, nspecies), flrwgt2(ld, md, nd, nspecies),  &
            flrwgt3(ld, md, nd, nspecies),   gewgt(ld, md, nd, nspecies),  &
               twgt(ld, md, nd, nspecies),    nwgt(ld, md, nd, nspecies),  &
           pfilter3(ld, md, nd, nspecies))

    allocate (pfilter(ld, md, nd), denom0(md, nd))

    allocate (dm_out(md, nd), time(md), dt(md))
    
    allocate (rkx(ld, md, nd), rkx2(ld, md, nd), rkperp2(ld, md, nd))
    allocate (rky(md, nd), rky2(md, nd), r0(md, nd))

    allocate (b_mag(ld),  bmaginv(ld), bgrad(ld),  gbd(ld), gbd0(ld), cvd(ld), &
               cvd0(ld), jacobian(ld),   gd2(ld), gd21(ld), gd22(ld), &
            gradrho(ld),        r(ld),    dr(ld)) 

    allocate (mmin(nd), mmax(nd))

  end subroutine allocate

  subroutine allocate_nez(ntz, ncntz)
    
    integer :: ntz, ncntz
!
! Allocate arrays that have a time index
!    
    allocate (timo(ntz), gamx(ntz), wrhx(ntz), wrhy(ntz), fluxe(ntz), qfluxe(ntz), &
             pcerr(ntz), wxsp(ntz), wakx(ntz), waky(ntz),  wenx(ntz),   wlux(ntz), &
              wdpx(ntz), wdvx(ntz),  upol(ntz),utor(ntz),phirms(ntz))

    allocate (phih0(ntz),uparh0(ntz),denh0(ntz),tparh0(ntz),tperph0(ntz), &
              qparh0(ntz),qperph0(ntz),veh0(ntz),uparph0(ntz))   

    allocate (wpfx(ntz, nspecies), fluxi(ntz, nspecies))

    allocate (grtmx(ntz, md, nd), wkif(ntz, md, nd),  wpif(ntz, md, nd), &
               wtif(ntz, md, nd),  psp(ntz, md, nd), phisq(ntz, md, nd))

    allocate (mgamx(md, nd, ntz))

    allocate (fluximn(ntz, md, nd, nspecies), qfluximn(ntz, md, nd, nspecies))
    if(epse > 0) allocate (fluxemn(ntz, md, nd), qfluxemn(ntz, md, nd))

    allocate (cx(ntz, 4*nd), cy(ntz, 4*md), cz(ntz, ld), czc(ntz, ld), czcn(ntz, ld), &
         czn(ntz, ld), cznn(ntz, ld))

    allocate (phi00(ntz, nd), den00(ntz, nd), upar00(ntz, nd), tpar00(ntz, nd), tperp00(ntz, nd))

    allocate (ct(ntz))
    allocate (drt(14, ntz))

    allocate (wenrk(ntz, md, nd), dketa(ntz, md, nd), dkpar(ntz, md, nd), &
              dktor(ntz, md, nd), dktdp(ntz, md, nd), wkups(ntz, md, nd), dktot(ntz, md, nd))
    
    allocate (timf(ncntz), utim(ncntz, md, nd))

  end subroutine allocate_nez

end module itg_data