@@ -581,10 +581,8 @@ contains
581581 -offset_z%beg:p + offset_z%end, nm), &
582582 intent(out) :: liutex_axis !< Liutex rigid rotation axis
583583
584- real(wp), dimension(nm) :: omega !< Vorticity
585-
586- character :: ivl = ' N' !< compute left eigenvectors
587- character :: ivr = ' V' !< compute right eigenvectors
584+ character, parameter :: ivl = ' N' !< compute left eigenvectors
585+ character, parameter :: ivr = ' V' !< compute right eigenvectors
588586 real(wp), dimension(nm, nm) :: vgt !< velocity gradient tensor
589587 real(wp), dimension(nm) :: lr, li !< real and imaginary parts of eigenvalues
590588 real(wp), dimension(nm, nm) :: vl, vr !< left and right eigenvectors
@@ -628,11 +626,6 @@ contains
628626 end do
629627 end do
630628
631- ! Compute vorticity
632- omega(1) = vgt(3, 2) - vgt(2, 3)
633- omega(2) = vgt(1, 3) - vgt(3, 1)
634- omega(3) = vgt(2, 1) - vgt(1, 2)
635-
636629 ! Call appropriate LAPACK routine based on precision
637630#ifdef MFC_SINGLE_PRECISION
638631 call cgeev(ivl, ivr, nm, vgt, nm, lr, li, vl, nm, vr, nm, work, lwork, info)
@@ -653,14 +646,16 @@ contains
653646 eigvec_mag = sqrt(eigvec(1)**2._wp &
654647 + eigvec(2)**2._wp &
655648 + eigvec(3)**2._wp)
656- if (eigvec_mag /= 0._wp ) then
649+ if (eigvec_mag .gt. sgm_eps ) then
657650 eigvec = eigvec/eigvec_mag
651+ else
652+ eigvec = 0._wp
658653 end if
659654
660655 ! Compute vorticity projected on the eigenvector
661- omega_proj = omega(1 )*eigvec(1) &
662- + omega(2 )*eigvec(2) &
663- + omega(3 )*eigvec(3)
656+ omega_proj = (vgt(3, 2) - vgt(2, 3) )*eigvec(1) &
657+ + (vgt(1, 3) - vgt(3, 1) )*eigvec(2) &
658+ + (vgt(2, 1) - vgt(1, 2) )*eigvec(3)
664659
665660 ! As eigenvector can have +/- signs, we can choose the sign
666661 ! so that omega_proj is positive
0 commit comments