subroutine cnstnts(old_qpar,old_qperp,tim,dt)
c
c Set up arrays and define values of constants.
c
      implicit none
      include 'itg.par'
      include 'itg.cmn'
      complex old_qpar(lz,mz,nz,nspecz),old_qperp(lz,mz,nz,nspecz)
      integer i,l,m,n,itime(7),lp,lm,idum,ntablex
      real tim,dt,z1,z2,z3
      character*30 x05abE
      external x05abE
      external x05aaE
c
      ldb=ld-1
c
c initialize arrays of m and n modes being used:
c
c nd		# of n modes
c md	# of m modes for each n
c
      malias=0
      nalias=0
      if(nffty .ne. 0 .and. nfftx .ne. 0) then
         malias=nffty
         nalias=nfftx
         if(md .eq. 0 .and. nd .eq. 0) then
            md=(nffty-1)/3+1
            nd=2*((nfftx-1)/3)+1
         endif
      endif

      if(nd.eq.-1) then
         nmin=0
         nmax=1
         nddm=1
         nddp=1
         nrr(1)=1
         nd=1
         goto 413
      endif
c
      if(mod(nd,2).eq.0) nd=nd+1
      if(nd.eq.1) then
         nmin=0
         nmax=0
         nddm=1
         nddp=1
         nrr(1)=0
      else
         nddm=nd/2
         nddp=nd/2+1
         nmin=-nd/2
         nmax=nd/2
         do n=1,nddp
            nrr(n)=n-1
         enddo
c     
         do n=1,nddm
            nrr(n+nddp)=n-nddm-1
         enddo
      endif
c
 413  continue
      if(ld .gt. lz) call aborter(6,
     &     'error: ld> lz (see itg.cmn and itg.par)')
      if(nd .gt. nz) call aborter(6,
     &     'error: nd> nz (see itg.cmn and itg.par)')
      if(md .gt. mz) call aborter(6,
     &     'error: md> mz (see itg.cmn and itg.par)')

      meq=-1	! set defaults in case there is no m=n=0 mode

      do n=1,nd
         if(md.gt.0) then
            mmax(n)=md-1
            mmin(n)=0
         else
            md=mlow-mhi+1
            mmin(n)=mlow
            mmax(n)=mhi
         endif
         do m=1,md
            mrr(m,n)=m+mmin(n)-1
c     
            if(mrr(m,n).eq.0.and.nrr(n).eq.0) then
               meq=m
               neq=n
            endif
c     
         enddo
      enddo

c Set up for real-to-complex y FFT's.
c md = # of independent ky modes (with ky>=0):
      mdspec=2*(md-1)           ! = # of ky points (without dealiasing)
      if(malias .eq. 0) then
         malias=4*(md-1)        ! overkill (more efficient to specify nffty)
      endif
      if(malias .gt. mffty) call aborter(6,
     &     'error: nffty or malias > mffty (see itg.in and itg.par)')
      
      madd=malias-mdspec
      rootm=sqrt(float(malias))
      
c     initialize trig tables for y FFT's:
      idum=0
      call csfftm(0,malias,1,1.0,idum,idum,idum,idum,tabley,idum,0)
c     
c     Set up for complex-to-complex x FFT's.
c     nd = # of independent kx modes (usually odd for symmetry around kx=0)
c     nalias = # of x points (with dealiasing):
      if(nd.eq.1) nalias=1
      if(nalias .eq. 0) then
         nalias=2*(nd-1)        ! overkill (more efficient to specify nfftx)
      endif
      if(nalias .gt. mfftx) call aborter(6,
     &     'error: nfftx or nalias > mfftx (see itg.in and itg.par)')
      
      nadd=nalias-nd
      rootn=sqrt(float(nalias))
      
c     initialize trig tables for x FFT's:
      ntablex=100+2*nzz
      call xmcfft(0,nalias,1,1.0,0.0,idum,idum,0.0,idum,idum,
     &     tablex,ntablex,1.0,idum)
c     
c     finish setting default m=0,n=0 mode if not already specified:
      if(meq .eq. -1) then
         neq=1
         meq=md+1
         if(meq .gt. mz) call aborter(6,'no room for m=0 mode')
      endif
c     
      pi=3.141592653589793
c     pi=abs(acos(-1.0))  ! machine independent form...

c     rkkin is the coefficient for the Hammett-Perkins
c     Landau damping model.
c     
c     rkkin = chi_1 * sqrt(2) * v_ti / c_s
c     = 2 * (2/pi*Ti/Te)**0.5,
c     
c     
      rkkin=2.*(2./pi)**0.5
      rkkperp=(2./pi)**0.5
      rkkperp2=(pi/(2.))**0.5
c     
c     zDq1 and zBq1 are the coefficients for the 4-pole Landau damping model
c     from Greg's PRL:
c     
c     zDq1 = D1 * sqrt(2) * Sqrt(TiovTe)
c     D1 = 2 sqrt(pi) / (3 pi -8)
c     zBq1 = (3 + 2*B1)*TiovTe
c     zEq1 = (3 + 2*B1)
c     B1 = (32 - 9 pi) / (6 pi - 16)
      
      zDq1 = 2.*sqrt(pi)/(3.0*pi-8.0)*sqrt(2.)
c     zBq1 = 3.+2.*(32.-9.*pi)/(6.*pi-16.)   bug in 4-moment closure
      zBq1 = (3.+2.*(32.-9.*pi)/(6.*pi-16.))
      zEq1 = 3.+2.*(32.-9.*pi)/(6.*pi-16.)
c     
      if(lin.eq.1) igradon=0
c     
      avgflag=0.0
      navgflag=0
      ncycle = 0
      ne=0
c     
c     take a radial point sqrt(L_s/L_n)*rho from the rational surface 
c     centered in the box and keep the time history of phi there.
c     
      if(shr .ne. 0.0) then
         lhistory=sqrt(1./abs(shr))*ld/x0
      else
         lhistory=ld/4
      endif
      if(iodd .eq. 2) lhistory=0.
      if(lhistory .gt. ld/4) lhistory=ld/4
c     lhistory=lhistory + ld/2+1  
c     
      lhistory=float(ld)*.6+1
c     changed mbeer (old way didn't work for large x0)
c     
c     when fortran does integer division, fractional remainders are
c     set to zero, i.e., 2/3 is 0, 5/4 is truncated to 1:
      nfreq=(nstp-1)/nez+1
      nfreqcnt=1
      if (nfreq.ne.0) then
         ntotal=nstp/nfreq
      else
         ntotal=1
      endif
c     
c     r(l) is the grid for the radial coordinate:
c     
      write(9,*) 'Theta_0:'
      do l=1,ld
         r(l)=2.*x0*float(l-1)/ldb-x0
      enddo
      do m=1,md
         do n=1,nd
            if (mrr(m,n).eq.0) then
               if (nmax.ne.0) then
                  r0(m,n)=z0*float(nrr(n))/nmax
               else
                  r0(m,n)=z0*float(nrr(n))
               endif
            else
               if (nmax.ne.0) then
                  r0(m,n)=z0*float(nrr(n))/nmax/float(mrr(m,n))
               else
                  r0(m,n)=z0*float(nrr(n))/float(mrr(m,n))
               endif
            endif
            write(9,*) m,n,' r0= ',r0(m,n)
         enddo
      enddo
c     
c     center of theta axis moved to zero, r0=theta0, mbeer
c     
c     Find cfl multipliers:
c     
c     cflx=float(ldb)/x0
c
c This line gives trouble if nd=1, but we don't run nonlinearly with nd=1,
c so just replace it with something similar:
c
c      cflx=shr*z0*float(nd)/(float(nd/2)*2*pi*y0)
      if(nd.eq.1) then
         cflx=shr*z0/(2.*pi*y0)
      else
         cflx=shr*z0*float(nd)/(float(nd/2)*2*pi*y0)
      endif
      cfly=float(mdspec)/(2.*pi*y0)
c     
c     bmag(l) is the magnitude of the magnetic field B normalized to B_0
c     bmaginv(l) is 1/B
c     bgrad is 1/B dB/dl    
c     
c     All in the circular concentric small aspect ratio limit for now.
c     
      do l=1,ld
         bmag(l)=1./(1.+eps*cos(r(l)))
         bmaginv(l)=1./bmag(l)
         bgrad(l)=epsn*eps/qsf*sin(r(l))*bmag(l)
      enddo   
c     
c     define ky grid
c     
      do n=1,nd
         do m=1,md
            rky(m,n)=(mrr(m,n)/y0) ! ky
            rky2(m,n)=(mrr(m,n)/y0)**2 ! ky**2
         enddo
      enddo
c     
c     define kx grid:
c     
c The current version (6.0.1.1) of Cray's fpp has a bug in the following
c nested loops, unless vectorization is turned off:
cfpp$ skip R
      do n=1,nd
         do m=1,md
c     
c     Need switch for m=0 modes.
c     
            do l=1,ld
               rkx(l,m,n)=(mrr(m,n)/y0)*((r(l)-r0(m,n))*shr
     &              -alpha*(sin(r(l))-sin(r0(m,n))))
               if (m.eq.meq) then
                  rkx(l,m,n)=mrr(m+1,n)/y0*(-r0(m+1,n))*shr
               endif
               if (md.eq.1) then
                  rkx(l,m,n)=1./y0
               endif
c     
c     local if qsf < 0
c     
               if (qsf.lt.0) rkx(l,m,n)=0.0
c     
               rkx(l,m,n)=rkx(l,m,n)*bmaginv(l)
               rkx2(l,m,n)=rkx(l,m,n)**2
            enddo
         enddo
      enddo
c     
cfpp$ noskip R
c     
c     omegad(l) is the w_d profile (for toroidal terms) without k_y
c
      do m=1,md
         do n=1,nd
            do l=1,ld
               omegad(l,m,n)=-2.*bmag(l)*epsn*rky(m,n)*( cos(r(l)) 
     &              +(shr*(r(l)-r0(m,n))
     &              -alpha*(sin(r(l))-sin(r0(m,n))))*sin(r(l)))
c     
c     special case for k_theta=0, k_r=const.     
c     
               if ((m.eq.meq).or.(md.eq.1)) then
                  omegad(l,m,n)=-2.*bmag(l)*epsn*rkx(l,m,n)*sin(r(l))
               endif
c     
c     local if qsf<0
c     
               if (qsf.lt.0) omegad(l,m,n)=-2.*epsn*rky(m,n)
            enddo
         enddo
      enddo
c     
      do l=2,ldb
         dr(l)=r(l+1)-r(l-1)
      enddo
      dr(1)=2.*(r(2)-r(1))
      dr(ld)=2.*(r(ld)-r(ldb))
c     
c     Twisted periodicity boundary conditions parameters:
c     
      if (iperiod.eq.2) then
         l_left=float(ld-1)/2.*(1.-xp/x0)+1.+.5
         l_right=float(ld-1)/2.*(1.+xp/x0)+1.+.5
         write(9,*) 'l_left, l_right: ', l_left, l_right
      endif
c     
      call flrinit
      if(nread.eq.1) then
      call wread(tim,dt)
      else
         do i=1,nspecies
            do n=1,nd
               do m=1,md
                  do l=1,ld
                     potential(l,m,n,i)=0.
                     density(l,m,n,i)=0.
                     u_par(l,m,n,i)=0.
                     t_par(l,m,n,i)=0.
                     q_par(l,m,n,i)=0.
                     t_perp(l,m,n,i)=0.
                     q_perp(l,m,n,i)=0.
c
c The next two lines prevent an uninitialized variable problem
c for a 3+1 model:
c
                     old_qpar(l,m,n,i)=0.
                     old_qperp(l,m,n,i)=0.
                     nl_density1(l,m,n,i)=0.
                     nl_density2(l,m,n,i)=0.
                     nl_upar1(l,m,n,i)=0.
                     nl_upar2(l,m,n,i)=0.
                     nl_tpar(l,m,n,i)=0.
                     nl_qpar(l,m,n,i)=0.
                     nl_tperp1(l,m,n,i)=0.
                     nl_tperp2(l,m,n,i)=0.
                     nl_tperp3(l,m,n,i)=0.
                     nl_qperp1(l,m,n,i)=0.
                     nl_qperp2(l,m,n,i)=0.
                  enddo
               enddo
            enddo
         enddo
      call pertb
      tim=0.0
      endif
c
c Get the time and date of this run:
c
      call x05aaE(itime)
      date=x05abE(itime)
c
      return
      end