subroutine itg(icontrol)
!
!     (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, 
!     and Gregory W. Hammett. ALL RIGHTS RESERVED.
!
! See itg.hlp for more info.
!
  use itg_data
  use in, only: input
  use flr_terms, only: flr, flrinit
  use nl, only: courant
  use step, only: timestep
  use diagnostics, only: prstep, volsq, posten, preenr
  implicit none

  complex, dimension(:,:,:,:), pointer, save :: old_n, old_upar, old_tpar, old_qpar, old_tperp, old_qperp       
  complex, dimension(:,:,:),   pointer, save :: old_pot, save_phi, old_e_density, old_e_p, old_e_r, old_e_t, phi_bat
  complex, dimension(:,:),     pointer, save :: omegaz

  real, dimension(:),   pointer, save :: hdt, k_tim

  real g_target,error,d_target
  real t_start
  integer m,n,ncyclea,ncycleb,icontrol(2)

  interface 
     subroutine baphi(phi)
       complex, dimension(:,:,:) :: phi
     end subroutine baphi
  end interface
  interface
     subroutine poisson(phi, a, b)
       complex, dimension(:,:,:) :: phi
       complex, dimension(:,:,:,:) :: a, b
     end subroutine poisson
  end interface
  interface
     subroutine kapav(a, b)
       complex, dimension(:,:,:) :: a, b
     end subroutine kapav
  end interface
  interface
     subroutine nonlin(a, b, c, d, e, f, phi)
       complex, dimension(:,:,:) ::  phi
       complex, dimension(:,:,:,:) :: a, b, c, d, e, f
     end subroutine nonlin
  end interface
  interface 
     subroutine nonline(a, b, c, d)
       complex, dimension(:,:,:) :: a, b, c, d
     end subroutine nonline
  end interface
  interface 
     subroutine filter(a, b, c, d, e, f, dt)
       complex, dimension(:,:,:,:) :: a, b, c, d, e, f
       real, dimension(:) :: dt
     end subroutine filter
  end interface

! input files:
!
! 12 RUNNAME.in     main namelist input file to control the run.
!  5 RUNNAME.resp   a copy of itg.res to continue from a previous run.
!
! output files:
!
!  7 RUNNAME.res    a complete dump of the common blocks.
!  9 RUNNAME.out    printout of various results
! 20 RUNNAME.fields phi(x,y) for particle code.
!
! **************************************************************
!
! Read namelists and set defaults
  call input

! Set up arrays, constants, and ic's
  call init
! initialize local arrays
  old_qpar    = 0.  ;  old_qperp   = 0.  ;  if(epse > 0.) old_e_density  = 0.
  icontrol(2)=1

! Allocate local arrays
  call allocate_local 
  
  if(ifilter == -1) t_start=time(md)
!
! ****************************************************
!                     Top of the loop
! ****************************************************
100 continue

  ncycle=ncycle+1
!
! Find L_T_crit:
!
  if(icrit /= 0 .and. ncycle >= 26) call critical(g_target,d_target,icontrol(2))
!
! The code advances the fields in time using a simple explicit 2-step
! Runge-Kutta algorithm that is second order accurate in time.
!
! Take the first step with hdt=0.5*dt
!
! Calculate FLR corrections to the potential
!
  call flr(potential)
  if (epse > 0.0) call baphi(potential)

! Evaluate non-linear terms if needed

  if(lin /= 1) then
     call nonlin(density, u_par, t_par, q_par, t_perp, q_perp, potential)
     if (epse > 0.0) call nonline(e_density, e_p, e_r, phi_ba)
  endif
!
! Simple-minded way to adjust the time step for nonlinear runs to the 
! radial Courant condition
!
  if(lin /= 1 .and. idt == 1) then
     if (.not.((nread == 1) .and. (ncycle < 3))) call courant
  endif
!     
  
  time = time + dt
  hdt =0.5*dt
!     
! Take a half time step.
!
  call timestep(old_n, old_upar, old_tpar, old_qpar, old_tperp, old_qperp, &
      density, u_par, t_par, q_par, t_perp, q_perp, old_e_density, &
      old_e_p, old_e_r, old_e_t, e_density, e_p, e_r, e_t, hdt)
!
! Solve Poisson equation
!
  if (epse > 0.0) call kapav(phi_ba, old_e_density)
  call poisson(old_pot, old_n, old_tperp)
!
!  If ifilter=1, include (hyper-)viscosity-like terms in each equation.
!   (n,v,tpar,qpar,T,q)
!     
  if (iabs(ifilter) == 1) then
     call filter(old_n, old_upar, old_tpar, old_qpar, old_tperp, old_qperp, hdt)     
  endif
!     
!*****************************************************************
!  Take the second step with a full dt to find values at (n+1)*dt
!*****************************************************************
!
! Evaluate various FLR corrections to phi. 
!
  call flr(old_pot)
  if (epse > 0.0) call baphi(old_pot)
!     
! Evaluate non-linear terms if needed.
! 
  if(lin /= 1) then
     call nonlin(old_n, old_upar, old_tpar, old_qpar, old_tperp, old_qperp, old_pot)
     if (epse > 0.0) call nonline(old_e_density, old_e_p, old_e_r, phi_ba)
  endif
!     
! Take a full time step.
!
  call timestep(density, u_par, t_par, q_par, t_perp, q_perp,&
       old_n,old_upar, old_tpar, old_qpar, old_tperp, old_qperp, e_density, &
       e_p, e_r, e_t, old_e_density, old_e_p, old_e_r, old_e_t, dt)
!
!     check particle conservation
!
!      if (mod(ncycle,nprnt) == 0) then 
!
!     save old phi_ba, so it's not overwritten
!
!         phi_bat = phi_ba

!         call baphi(potential)
!         call pcons
!
!     reload old phi_ba
!
!         phi_ba = phi_bat

!      endif
!
! Solve Poisson equation: 
!
  if (epse > 0.0) call kapav(phi_ba,e_density)
  call poisson(potential,density,t_perp)
!
! If ifilter = 1, include (hyper-)viscosity-like terms in each equation.
!   (n,v,tpar,qpar,T,q)
!     
  if (iabs(ifilter) == 1) then
     call filter(density, u_par, t_par, q_par, t_perp, q_perp, dt)
  endif
!
! Keep time history of the potential.
!
  if(mod(ncycle, nfreq) == 0) then
     nfreqcnt=nfreqcnt+1
     utim(nfreqcnt,:,:)=potential(lhistory,:,:)
     timf(nfreqcnt)=time(md)
  endif
!
! Check the energy balance.
!
  ncyclea=ncycle+1
  ncycleb=ncycle-1
 
  if(mod(ncyclea,nprnt) == 0) call preenr  
  if(ncycle /= 1 .and. mod(ncycleb,nprnt) == 0) call posten

!     
!     Save some results for printing later.
!     
  if(mod(ncycle,nprnt) == 0 .or. (ncycle == 1)) then
     ne = ne + 1
     call prstep(t_start)
     if(icrit /= 0 .and. ncycle >= 26) call tc_out(g_target,d_target,error,icontrol)
!
! if ifilter = -1, use time-dependent viscosity-like terms in each equation.
!   (n,v,tpar,qpar,T,q) proportional to sqrt(int(dt chi))
!      
     if(ifilter == -1) call rdiffset
     
  endif
  if(mod(ncycle+1, nprnt) == 0 .and. icrit /= 0) save_phi=potential
     
  if(ntrace == 0) goto 1000
  if(mod(ncycle,nprnt) == 0) then
     if(ntrace == 1) then
        if (epse > 0.0) then
           write(*,2000) ncycle,dt(md),wpfx(ne,1),qfluxe(ne), &
                wpfx(ne,1)/qfluxe(ne),fluxe(ne),fluxi(ne,1)
        else   
           if(ifilter == -1) then
              write(*,2002) ncycle,rdiff(1,1),-wpfx(ne,1)/eta(1)
              write(123,*) timo(ne),rdiff(1,1),-wpfx(ne,1)/eta(1), &
                   chi_int,1.5*rdiff(1,1)-wpfx(ne,1)/eta(1)
           else
              write(*,2001) ncycle,dt(md),wpfx(ne,1)
           endif
        endif
     else
        if (epse > 0.0) then
           write(59,2000) ncycle,dt(md),wpfx(ne,1),qfluxe(ne), &
                wpfx(ne,1)/qfluxe(ne),fluxe(ne),fluxi(ne,1)
        else   
           write(59,2001) ncycle,dt(md),wpfx(ne,1)
        endif
     endif
  endif
2000 format(' n: ',i4,' dt: ',f6.4,' Q_i: ',1pe10.2,' Q_e: ',1pe10.2, &
          ' Q_i/Q_e: ',1pe10.2, ' Gamma_e: ',1pe10.2,' Gamma_i: ',1pe10.2)
2001 format(' n: ',i4,' dt: ',f6.4,' Q_i: ',1pe10.2,' adiabatic electrons')
2002 format(' n: ',i4,' D: ',1pe10.2,' chi: ',1pe10.2)
1000 continue
      
  if(icontrol(1) == 1) goto 200
  if(ncycle < nstp) goto 100
200 continue
  
!
!     take out B variation of kperp rho for post gamma/k^2 calculation
!      
  do n=1,nd
     do m=1,md
        rkperp2(:,m,n)=rkperp2(:,m,n)/bmaginv(:)**2
     enddo
  enddo
  
contains

  subroutine tc_out(g_target,d_target,err,icontrol)

!
! Write critical temperature gradient info
! or D_R info
!
    complex, dimension(:,:), pointer, save :: omega0
    real,    dimension(:,:), pointer, save :: e1, g1, dml, e_crit, D_old
    real,    dimension(:,:), allocatable :: e2, g2
    complex w0,w1
    complex phi1,phi0
    real err
    real g_rate
    real D_max,g_max
    real g_target,d_target
    real dmix
    integer :: done = 0
    integer i,l,m,n,init,icontrol(2),lsave,m_min
    
    if(.not.associated(omegaz)) allocate (omegaz(md, nd), e1(md, nd), D_old(md, nd), &
         g1(md, nd), dml(md, nd), e_crit(md, nd), k_tim(md), omega0(md, nd))
    m_min=2
    if(md == 1) m_min=1
    
    if(icontrol(2) == 2) init=1
    if(done == 3) done=0
    icontrol(2)=3
    
    lsave=ldb/2+2
    if(init == 1) then
       init=0
       do m=m_min,md
          D_old(m,:)=0.
          omegaz(m,:)=1.
          dml(m,:)=1.
       enddo
    else
       g_max=-1.e10
       do n=1,nd
          do m=m_min,md
             phi0=save_phi(lsave,m,n)
             phi1=potential(lsave,m,n)
             omega0(m,n)=omegaz(m,n)
             omegaz(m,n)=ii*log(phi1/phi0)/dt(m)
             g_max=max(g_max,aimag(omegaz(m,n)))
             D_old(m,n)=dml(m,n)
          enddo
       enddo
       
       call dm_calc(potential,omegaz,dml,D_max)
    endif
    
    if(icrit == 1) then
       err=0.
       do n=1,nd
          do m=m_min,md
             g_rate=aimag(omegaz(m,n))
             w0=omega0(m,n)
             w1=omegaz(m,n)
             err=max(err,sqrt(cabs((w1-w0)/(w1+w0))))
             err=max(err,abs((g_rate-g_target)/(g_target+g_rate)))
          enddo
       enddo
    endif
    if(icrit == 2) then
       err=0.
       do n=1,nd
          do m=m_min,md
             dmix=dml(m,n)
             err=max(err,abs((dmix-d_target)/(d_target+dmix)))
             err=max(err,abs((dmix-D_old(m,n))/(dmix+D_old(m,n))))
          enddo
       enddo
    endif
    if(err < tol .or. ncycle == nstp) done=done+1
    
    if(done == 1) then
       done=2
       if(icontrol(1) == 0) open(unit=1,file=runname(1:lrunname)//'.dmix')
!       if(icrit == 1) write(*,*) D_max,g_max,etak(m_min,1,1), &
!            ncycle,err,g_target,d_target
       if(icrit == 2) then
!           write(*,*) D_max,ncycle,err,d_target
!           write(*,*) (rdiff(j),j=1,n_ktheta)
          done=3
       endif
       g_target=gamma_0/2.
       do m=m_min,md
          g1(m,:)=aimag(omegaz(m,:))
          e1(m,:)=etak(m,:,1)
       enddo
       if(ncycle == nstp) done=3
    endif
    
    if(done == 3) then
       if(icrit == 1) then
!            write(*,*) D_max,g_max,etak(m_min,1,1),ncycle,err,g_target
          allocate (e2(md, nd), g2(md, nd))
          do m=m_min,md
             g2(m,:)=aimag(omegaz(m,:))
             e2(m,:)=etak(m,:,1)
             where (g2(m,:)-g1(m,:) /= 0) 
                e_crit(m,:)=e2(m,:)-g2(m,:)*(e2(m,:)-e1(m,:))/(g2(m,:)-g1(m,:))
             elsewhere
                e_crit(m,:)=-1.e10
             end where
          enddo
          deallocate (e2, g2)
          write(*,*)
          do n=1,nd
             do m=m_min,md
                if(e_crit(m,n) == -1.e10) then
                   write(*,2002) rky(m,n),r0(m,n)
                else
                   write(*,2001) time(m),rky(m,n),r0(m,n),&
                        e_crit(m,n),e_crit(m,n)/epsn
                endif
             enddo
          enddo
          write(*,*) 
          icontrol(1)=1
       endif
       if(icrit == 2) then
          do m=m_min,md               
             dm_out(m,:)=dml(m,:)+rdiff(m,:)
          enddo
          write(*,*) 
          do n=1,nd
             do m=m_min,md
                write(*,2000) time(m),rky(m,n),r0(m,n),dm_out(m,n)
             enddo
             write(*,*) 
          enddo
          icontrol(1)=1
       endif
       return
    endif
    
2000 format(' time= ',f7.2,'  k_theta= ',f6.4,'  theta_0= ',f8.6,'  D_ML= ',f8.5)
2001 format(' time= ',f7.2,'  k_theta= ',f6.4,'  theta_0= ',f6.4,'  Ln/LT_crit= ',f8.5,' R/LT_crit= ',f8.5)
2002 format('Failed to find LT_crit for k_theta= ',f6.4,'  theta_0= ',f6.4)
    
  end subroutine tc_out
  
  subroutine dm_calc(a,omegaz,dml,D_max)
  
    complex, dimension(:,:,:) ::  a
    complex, dimension(:,:) :: omegaz
    real, dimension(:,:) :: dml
    real :: norm, D_max, avekp2
    integer :: l, m, n
    
    D_max=0.
    do n=1,nd
       do m=1,md
          avekp2=0.0
          norm=0.0
          do l=1,ldb
             avekp2=avekp2+rkperp2(l,m,n)*cabs(a(l,m,n))**2
             norm=norm+cabs(a(l,m,n))**2
          enddo
          if (norm /= 0.0 .and. avekp2 /= 0.) then
             avekp2=avekp2/norm
             dml(m,n)=aimag(omegaz(m,n))/avekp2
          else
             dml(m,n)=0.
          endif
          D_max=max(D_max,dml(m,n))
       enddo
    enddo
  
  end subroutine dm_calc

  subroutine critical(g_target,d_target,icontrol2)
!
! Find the critical temperature gradient or the diffusion needed 
! to stabilize the mode.
!
! gamma_0 should be made a function of k_theta.
!
    real, dimension(:,:), pointer, save :: p2, etak_0, d_0, k_p2

    real g_target, d_target
    real dum, d_fac, e_fac, avekp2

    integer i, l, m, n, init, icontrol2, m_min
    save init

    if(.not.associated(p2)) allocate (p2(md, nd), etak_0(md, nd), d_0(md, nd), k_p2(md, nd))
!
! Assume m=1 is ky=0 mode unless md=1; don't calculate anything for 
! ky=0 mode here:
!
    m_min=2
    if(md == 1) m_min=1
    
    if(icontrol2 == 1) then
       init=icontrol2
       icontrol2=2
    endif
    if(init == 1) then
       init=2
       g_target=gamma_0
       d_target=rdiff_0
       call volsq(potential,potential,dum,p2)
       do m=m_min,md
          etak_0(m,:)=etak(m,:,1)
          d_0(m,:)=rdiff(m,:)
          k_tim(m)=time(m)
          k_p2(m,:)=p2(m,:)
       enddo
    else if(init == 2) then
       call volsq(potential,potential,dum,p2)
       if(icrit == 1) then
          if(gamma_0 /= g_target) then
             do m=m_min,md
                etak_0(m,:)=etak(m,:,1)
                k_tim(m)=time(m)
                k_p2(m,:)=p2(m,:)
             enddo
             gamma_0=g_target
          endif
          do m=m_min,md
             e_fac=exp(2.*gamma_0*(time(m)-k_tim(m)))
             etak(m,:,1)=etak_0(m,:)/sqrt(p2(m,:)/k_p2(m,:)/e_fac)
             etak_par(m,:,1)=etak(m,:,1)
          enddo
          if(nspecies > 1) then
             do i=2,nspecies
                do m=m_min,md
                   etak(m,:,i)=etak(m,:,1)
                   etak_par(m,:,i)=etak_par(m,:,1)
                enddo
             enddo
          endif
          call flrinit(0)
       else if(icrit == 2) then
          if(rdiff_0 /= d_target) then
             do m=m_min,md
                d_0(m,:)=rdiff(m,:)
                k_tim(m)=time(m)
                k_p2(m,:)=p2(m,:)
             enddo
             rdiff_0=d_target
          endif
          do n=1,nd
             do m=m_min,md
                avekp2=0.
                dum=0.
                do l=1,ldb
                   avekp2=avekp2+cabs(potential(l,m,n))**2*rkperp2(l,m,n)
                   dum=dum+cabs(potential(l,m,n))**2
                enddo
                if(dum /= 0.) then
                   avekp2=avekp2/dum
                else
                   avekp2=0.
                endif
                d_fac=exp(2.*rdiff_0*avekp2*(time(m)-k_tim(m)))
                rdiff(m,n)=d_0(m,n)*sqrt(p2(m,n)/k_p2(m,n)/d_fac)
             enddo
          enddo
       endif
    else 
       return      
    endif
    
  end subroutine critical

  subroutine rdiffset
!      
! This routine calculates the diffusion associated with particle noise
! in GK PIC simulations.  The code evaluates an analytical derivation 
! by Greg Hammett.  The parameters ifilter, N_p, n_zpic, and bmax are 
! described in the input file (itg.in).  
!      
      
    integer m, n
    real tmp, fa_par, fac
  
    fa_par=2*pi*qsf/epsn/float(n_zpic)   ! a_par/L_n
  
    do n=1,nd
       do m=1,md
          tmp=0.5*sqrt(pi/2)*eta(1)**2/float(n_p)*fa_par*chi_int*bmax/(2+bmax)
!     first time step hack
          if(rdiff(m,n) == 0.) then
             rdiff(m,n)=tmp
          else
             fac=1+sqrt(2/pi)/fa_par/(rdiff(m,n)-wpfx(ne,1)/eta(1)/1.5) &
                  *2/bmax*(2+bmax)/(1+bmax)
             rdiff(m,n)=tmp*alog(fac)
          endif
       enddo
    enddo
    
  end subroutine rdiffset

  subroutine allocate_local

    allocate (hdt(md))

    allocate (old_n(ld,md,nd,nspecies), old_upar(ld,md,nd,nspecies), &
         old_tpar(ld,md,nd,nspecies), old_qpar(ld,md,nd,nspecies), &
         old_tperp(ld,md,nd,nspecies), old_qperp(ld,md,nd,nspecies))
    
    allocate (old_pot(ld,md,nd), save_phi(ld,md,nd))

    if(epse > 0.) &
         allocate (old_e_density(kdpass, md, nd), old_e_p(kdpass, md, nd), &
         old_e_r(kdpass, md, nd), old_e_t(kdpass, md, nd), phi_bat(kdpass, md, nd))
    
  end subroutine allocate_local

end subroutine itg