th_bishop Subroutine

  pure subroutine th_bishop(rpgrad, th_bish, nth)
    use constants, only: twopi
    implicit none
    integer, intent (in) :: nth
    real, dimension (-ntgrid:, :), intent (in) :: rpgrad
    real, dimension (-ntgrid:), intent (out) :: th_bish
    real, dimension (-nth:nth) :: magrp
    real, dimension (-nth:nth, 2) :: tvec
    integer :: i
    real :: tvec_dot_gradR

    magrp = hypot(rpgrad(-nth:nth, 1), rpgrad(-nth:nth, 2))

    ! tvec is the tangent vector in the surface of the poloidal plane
    tvec(:, 1) = rpgrad(-nth:nth, 2) / magrp
    tvec(:, 2) = -rpgrad(-nth:nth, 1) / magrp

    ! need to find grad(R) and dot it with tvec.  but grad(R) is
    ! simple -- it is just (1,0) so simply return acos(tvec(:,1)) with
    ! corrections for left-hand quadrants
    do i = -nth, nth
       tvec_dot_gradR = min(1., max(-1., tvec(i, 1)))
       if (tvec(i, 2) > 0.) then
          th_bish(i) = twopi - acos(tvec_dot_gradR)
       else
          th_bish(i) = acos(tvec_dot_gradR)
       end if
    end do
  end subroutine th_bishop