subroutine tc_out(save_phi,g_target,d_target,
     .     err,time,dt,icontrol)
c
c Write critical temperature gradient info
c or D_R info
c
      implicit none
      include 'itg.par'
      include 'itg.cmn'

      complex save_phi(lz,mz,nz),w0,w1
      complex omegaz(mz,nz),phi1,phi0,omega0(mz,nz)
      real aveth,norm,err,e2(mz,nz),g2(mz,nz),e_crit(mz,nz)
      real e1(mz,nz),g1(mz,nz),dummy(mz,nz),d_tot,k_avg,g_rate
      real dml(mz,nz),D_max,g_max,dk_theta
      real g_target,d_target,D_max_old(2),dt(mz),time(mz)
      real D_old(mz,nz),dmix
      integer i,l,m,n,done,init,icontrol(2),lsave,m_min
      namelist/crit_stuff/D_max,g_max,shr,qsf,
     .     eps,alpha,epsn,etak,TiovTe,rmass,charge,n_I,
     .     nu_ii,rky,omegaz,dml,ncycle,rho
      namelist/gradti/e_crit
      data done/0/
      save D_max_old,done,e1,g1,dml,omegaz

      m_min=2
      if(md.eq.1) m_min=1

      if(icontrol(2).eq.2) init=1
      if(done.eq.3) done=0
      icontrol(2)=3

      lsave=ldb/2+2
      if(init.eq.1) then
         init=0
         D_max_old(2)=2.e10
         D_max_old(1)=1.e10
         do n=1,nd
            do m=m_min,md
               D_old(m,n)=0.
               omegaz(m,n)=1.
               dml(m,n)=1.
            enddo
         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.eq.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.eq.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.lt.tol .or. ncycle.eq.nstp) done=done+1

      if(done.eq.1) then
         done=2
         if(icontrol(1).eq.0) 
     >        open(unit=1,file=runname(1:lrunname)//'.dmix')
         write(1,crit_stuff)
c         if(icrit.eq.1) write(*,*) D_max,g_max,etak(m_min,1,1),
c     .        ncycle,err,g_target,d_target
         if(icrit.eq.2) then
c           write(*,*) D_max,ncycle,err,d_target
c           write(*,*) (rdiff(j),j=1,n_ktheta)
            done=3
         endif
         g_target=gamma_0/2.
         do n=1,nd
            do m=m_min,md
               g1(m,n)=aimag(omegaz(m,n))
               e1(m,n)=etak(m,n,1)
            enddo
         enddo
         if(ncycle.eq.nstp) done=3
      endif

      if(done.eq.3) then
         write(1,crit_stuff)
         if(icrit.eq.1) then
c            write(*,*) D_max,g_max,etak(m_min,1,1),ncycle,err,g_target
            do n=1,nd
               do m=m_min,md
                  g2(m,n)=aimag(omegaz(m,n))
                  e2(m,n)=etak(m,n,1)
                  if(g2(m,n)-g1(m,n).ne.0) then
                     e_crit(m,n)=e2(m,n)-g2(m,n)*
     .                    (e2(m,n)-e1(m,n))/(g2(m,n)-g1(m,n))
                  else
                     e_crit(m,n)=-1.e10
                  endif
               enddo
            enddo
            write(*,*)
            do n=1,nd
               do m=m_min,md
                  if(e_crit(m,n).eq.-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(*,*) 
            write(1,gradti)
            icontrol(1)=1
         endif
         if(icrit.eq.2) then
            do n=1,nd
               do m=m_min,md               
                  dm_out(m,n)=dml(m,n)+rdiff(m,n)
               enddo
            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)
      return
      end