anchor.diff:                  deallocate(phinew(i)%p)
                   allocate(phinew(i)%p(ry(i)*ry(i)*ra(i)*2))
                end if
  -             call dphi_right(ry(i), n(i), ry(i+1), ry(i), n(i), ry(i+1), ra(i), ra(i+1), phinew(i+1)%p, crA(pa(i)), crnew(i)%p, crnew(i)%p, phinew(i)%p)
  +             call dphi_right(ry(i), n(i), ry(i+1), ry(i), n(i), ry(i+1), ra(i), &
  +             ra(i+1), phinew(i+1)%p, crA(pa(i)), crnew(i)%p, crnew(i)%p, phinew(i)%p)
             end if
          end if
          if ( dir > 0 ) then
  
positive.code:              !In this case, we have to put S(i+1) backwards in time (heh)
               call zqr(n(i)*ry(i), ry(i+1), curcr, R) 
  
               call zphi_left(ry(i), n(i), ry(i+1), ry(i), n(i), ry(i+1), ra(i), ra(i+1), phinew(i)%p, crA(pa(i)), curcr, curcr, phitmp)
               call zinit_sfun(ry(i+1), ry(i+1), ra(i+1), ry(i+1), ry(i+1), phitmp, phinew(i+1)%p)
               !anorm = znormest(ry(i+1)*ry(i+1), 4, zsfun_matvec, zsfun_matvec_transp)
               anorm = 1d0
  
negative.code:           if ( i < d ) then
               call zqr(ry(i)*n(i),ry(i+1),crnew(i)%p,R)
               !The size of the updated s would be ry(i+1) -> 
               call zphi_left(ry(i), n(i), ry(i+1), ry(i), n(i), ry(i+1), ra(i), &
               ra(i+1), phinew(i)%p, crA(pa(i)), crnew(i)%p, crnew(i)%p, phitmp)
               call zinit_sfun(ry(i+1), ry(i+1), ra(i+1), ry(i+1), ry(i+1), phitmp, phinew(i+1)%p)
               !anorm = znormest(ry(i+1)*ry(i+1),4,zsfun_matvec,zsfun_matvec_transp)
               anorm = 1d0
  