subroutine pertb
c
c This subroutine inserts the perturbation.
c
      implicit none
      real a,b,ax,g05caE,c1,c2,c3,c4,c5,c6,p00,in
      integer i,l,m,n
      include 'itg.par'
      include 'itg.cmn'
c
      write(9,302)  ncycle
302   format(i7,'*****  perturbation inserted ***************')
c
      do i=1,nspecies
         do n=1,nd
            do m=1,md
               do l=1,ld
                  potential(l,m,n)=0.
                  density(l,m,n,i)=0.
                  u_par(l,m,n,i)=0.
                  t_par(l,m,n,i)=0.
                  t_perp(l,m,n,i)=0.
               enddo
            enddo
         enddo
      enddo
      do m=1,md
         do n=1,nd
            phiav_old(m,n)=0.0
         enddo
      enddo   
c
c Choose parity of run.
c
      a=1.
      b=1.
      if(iodd.eq.1) a=0.	!iodd = 1 -> odd modes only
      if(iodd.eq.2) b=0.	!iodd = 2 -> even modes only
c     
      do n=1,nd
         do m=1,md
            do l=2,ldb
               p00=pmag*(b*r(l)+a)*exp(-r(l)**2)
               potential(l,m,n)=cmplx(p00,-p00)
            enddo
         enddo
      enddo
c
c  randomize: don't want to multitask this section
c

CMIC$ GUARD

      call g05cbE(iseed)        !random ic's, initialize seed
      
      do n=1,nd
         do m=1,md
            a=g05caE(a)
            b=g05caE(b)
            do l=2,ldb
               potential(l,m,n)=
     .              cmplx(a*real(potential(l,m,n)),
     .              b*aimag(potential(l,m,n)))
            enddo
         enddo
      enddo
CMIC$ END GUARD

c     
c     gaussian blob in center of box
c     
c       a=pi*y0/4.*pmag
c       do l=1,ldb
c          do n=1,nd
c             us(l,1,n,1)=0.0
c             uc(l,1,n,1)=a*exp(-(a**2/2.*rkx2(ldb/2,1,n)))
c             do m=2,md
c               us(l,m,n,1)=0.0
c               uc(l,m,n,1)=2.*a*exp(-(a**2/2.*(rkx2(ldb/2,m,n)
c      &                     +rky2(m,n))))
c             enddo  
c          enddo   
c       enddo   

c
c the following code for initializing the density needs to be upgraded
c to use the new FLR approximations:
c
c     
c             bill suggests not doing this
      do i=1,nspecies
         do n=1,nd
            do m=1,md
               do l=1,ldb
                  density(l,m,n,i)=potential(l,m,n)
               enddo
            enddo
         enddo
      enddo

c iphi00=9  = set phi=0 for all Fourier components, not just (m=0,n=0).
c (no electric field, for neoclassical comparisons):
c Also, for this neoclassical test case, initialize a cos(k_r*r) 
c perturbation, which is constant along the field line, either in
c the density (pmag > 0) or the temperatures (pmag < 0).
c
      if((iphi00 .eq. 9).or.((md.eq.1).and.(mrr(1,1).eq.0))) then
         do i=1,nspecies
            do n=1,nd
               do m=1,md
                  do l=1,ld
                     density(l,m,n,i)=0.
                     potential(l,m,n)=0.
                  enddo
               enddo
            enddo
         enddo

         if(pmag .gt. 0.0) then
            do i=1,nspecies
               do n=1,nd
                  do m=1,md
                     do l=1,ld
c bumpy equilibria
      c1=1.0
      c2=0.
      c3=8*c2-2*c1
      c4=0.
      c5=0.
      c6=0.
c      density(l,m,n,i)=rkkperp2*eps*( (c1-2*c2)*sin(r(l))-eps/2*(c1-c2)*
c     &      sin(2*r(l)) )+eps*(c5-c4)*cos(r(l))+c6
c      u_par(l,m,n,i)=c1*b_mag(l)
c      q_perp(l,m,n,i)=-c1*b_mag(l)+c2*b_mag(l)**2
c      q_par(l,m,n,i)=c3*b_mag(l)-2*c2*b_mag(l)**2
c      t_par(l,m,n,i)=-rkkperp2*eps*( (c1-2*c2)*sin(r(l))-eps/2*(c1-c2)*
c     &      sin(2*r(l)))  + c4
c      t_perp(l,m,n,i)=-rkkperp2*eps*( (c1-2*c2)*sin(r(l))-eps/2*(c1-3*c2)*
c     &      sin(2*r(l)))-c4*eps*(cos(r(l))+eps)+c5*bmaginv(l)
c      potential(l,m,n)=0.0
c
c     good flow initialization:
      density(l,m,n,i)=pmag*rkx(l,m,n)*(1.-ii)
      potential(l,m,n)=pmag/rkx(l,m,n)*(1.-ii)
      u_par(l,m,n,i)=(pmag*(1-ii)
     &     -2*ii*rkx(l,m,n)*qsf*cos(r(l))*potential(l,m,n))*b_mag(l)
                     enddo
                  enddo
               enddo
            enddo
         else
            do i=1,nspecies
               do n=1,nd
                  do m=1,md
                     do l=1,ld
                        t_perp(l,m,n,i)=abs(pmag)
                        t_par(l,m,n,i)=abs(pmag)
                     enddo
                  enddo
               enddo
            enddo
         endif
      endif
c     
c           
c Enforce the reality condition for the ky=0 components:
c
      if(nd.ne.1) then
         do i=1,nspecies
            do n=1,nddm
               do l=1,ldb
                  potential(l,1,nddp+n)=conjg(potential(l,1,nddp+1-n))
                  density(l,1,nddp+n,i)=conjg(density(l,1,nddp+1-n,i))
                  u_par(l,1,nddp+n,i)=conjg(u_par(l,1,nddp+1-n,i))
                  t_par(l,1,nddp+n,i)=conjg(t_par(l,1,nddp+1-n,i))
                  t_perp(l,1,nddp+n,i)=conjg(t_perp(l,1,nddp+1-n,i))
               enddo
            enddo
         enddo
c     
c     This could be set to any real number to satisfy the 
c     reality condition.  
c     
         do i=1,nspecies
            do l=2,ldb
               potential(l,1,1)=0.
               density(l,1,1,i)=0.
               u_par(l,1,1,i)=0.
               t_par(l,1,1,i)=0.
               t_perp(l,1,1,i)=0.
            enddo
         enddo
      endif

      return
      end