subroutine critical(g_target,d_target,time,icontrol2)
c
c Find the critical temperature gradient or the diffusion needed 
c to stabilize the mode.
c
c gamma_0 should be made a function of k_theta.
c
      implicit none
      include 'itg.par'
      include 'itg.cmn'

      real p2(mz,nz),g_target,d_target,time(mz)
      real aveth,dum,d_fac,e_fac,avekp2
      real rkp2,etak_0(mz,nz),d_0(mz,nz)
      real k_tim(mz),k_p2(mz,nz)
      integer i,l,m,n,init,icontrol2,m_min
      save init,k_p2,k_tim,etak_0,d_0
c
c Assume m=1 is ky=0 mode unless md=1; don't calculate anything for 
c ky=0 mode here:
c
      m_min=2
      if(md.eq.1) m_min=1

      if(icontrol2.eq.1) then
         init=icontrol2
         icontrol2=2
      endif
      if(init.eq.1) then
         init=2
         g_target=gamma_0
         d_target=rdiff_0
         call volsq(potential,potential,dum,p2)
         do n=1,nd
            do m=m_min,md
               etak_0(m,n)=etak(m,n,1)
               d_0(m,n)=rdiff(m,n)
               k_tim(m)=time(m)
               k_p2(m,n)=p2(m,n)
            enddo
         enddo
      else if(init.eq.2) then
         call volsq(potential,potential,dum,p2)
         if(icrit.eq.1) then
            if(gamma_0.ne.g_target) then
               do n=1,nd
                  do m=m_min,md
                     etak_0(m,n)=etak(m,n,1)
                     k_tim(m)=time(m)
                     k_p2(m,n)=p2(m,n)
                  enddo
               enddo
               gamma_0=g_target
            endif
            do n=1,nd
               do m=m_min,md
                  e_fac=exp(2.*gamma_0*(time(m)-k_tim(m)))
                  etak(m,n,1)=etak_0(m,n)/sqrt(p2(m,n)/k_p2(m,n)/e_fac)
                  etak_par(m,n,1)=etak(m,n,1)
               enddo
            enddo
            if(nspecies.gt.1) then
               do i=2,nspecies
                  do n=1,nd
                     do m=m_min,md
                        etak(m,n,i)=etak(m,n,1)
                        etak_par(m,n,i)=etak_par(m,n,1)
                     enddo
                  enddo
               enddo
            endif
            call flrinit(0)
         else if(icrit.eq.2) then
            if(rdiff_0.ne.d_target) then
               do n=1,nd
                  do m=m_min,md
                     d_0(m,n)=rdiff(m,n)
                     k_tim(m)=time(m)
                     k_p2(m,n)=p2(m,n)
                  enddo
               enddo
               rdiff_0=d_target
            endif
            do n=1,nd
               do m=m_min,md
                  aveth=0.
                  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.ne.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
         
      return
      end