anchor.diff:      call Snsourceset(    numcells,a,x,dx,constq,qtype,qr,ql,sigt,siga,SnBCs,&
                            psiBCl,psiBCr,numangs)!Set qs & psiBCs
       iter=0
  -    do while (error>tol)
  +    flitermore = .true.
  +    do while (flitermore)
         iter=iter+1
         call FESn(         numcells,phiSnl,phiSnr,psil,psir,x,dx,qr,ql,sigt,siga,mu,&
                            psiBCl,psiBCr,numangs ) !Solve an Sn sweep
         call FESncalcphi(  numcells,phiSnl,phiSnr,psil,psir,wgts,numangs ) !Calc phi from psi
         if(iter>1) call FESnerror( numcells,phiSnl,phiSnr,phiSnlold,phiSnrold,error ) !Check converge?
         call FESnNewToOld( numcells,phiSnl,phiSnr,phiSnlold,phiSnrold )
  +      if(error<tol) flitermore = .false.
       enddo
   !    print *,"iteration: ",iter
   !    call cpu_time(toc)
  
positive.code: 
        if(iter>1) call FESnerror( numcells,phiDSAl,phiDSAr,phiDSAlold,phiDSArold,error ) !Converge?
        call FESnNewToOld( numcells,phiDSAl,phiDSAr,phiDSAlold,phiDSArold )
      enddo
  !    print *,"iteration: ",iter
  !    call cpu_time(toc)
  
negative.code:     !init for diff
      call fluxinitcont(   numcells,phiDSAdiff,phistart ) !create phi=phistart
      iter=0
      flitermore = .true.
      do while (flitermore)
        iter=iter+1
        !Step 1, Sn solve for half step solution of flux
        call FESn(         numcells,phiDSAlold,phiDSArold,psil,psir,x,dx,qr,ql,sigt,siga,mu,&
  