gradient Subroutine

  subroutine gradient(self, rgrid, theta, grad, char, rp, nth_used, ntm, use_bishop)
    use splines, only: inter_d_cspl
    implicit none
    class(abstract_geo_type), intent(in) :: self
    integer, intent (in) :: nth_used, ntm
    character(1), intent (in) :: char
    real, dimension (-ntm:), intent (in)  :: rgrid, theta
    real, dimension (-ntm:,:), intent (out) :: grad
    real, intent(in) :: rp
    logical, intent(in) :: use_bishop
    real, dimension(:, :, :), allocatable :: darr
    real, dimension (1) :: p_eqpsi, dp_deqpsi
    integer :: i
    if (use_bishop) then
       select case(char)
       case('B')
          darr = self%dbbish
       case('P', 'R')
          darr = self%dpbish
       case('T')
          darr = self%dtbish
       end select
    else
       select case(char)
       case('B')
          darr = self%dbcart
       case('P', 'R')
          darr = self%dpcart
       case('T')
          darr = self%dtcart
       end select
    end if

    do i = -nth_used, nth_used
       grad(i, 1) = self%eqitem(rgrid(i), theta(i), darr(:,:,1), 'R')
       grad(i, 2) = self%eqitem(rgrid(i), theta(i), darr(:,:,2), 'Z')
    end do

    if(char == 'T' .and. .not. self%has_full_theta_range) then
       where (theta(-nth_used:nth_used) < 0.0)
          grad(-nth_used:nth_used, 1) = -grad(-nth_used:nth_used, 1)
          grad(-nth_used:nth_used, 2) = -grad(-nth_used:nth_used, 2)
       end where
    end if

    !     to get grad(pressure), multiply grad(psi) by dpressure/dpsi
    ! Could use self%dpfun if took pbar rather than rp?
    if (char == 'R') then !Only actually used if bishop == 0
       call inter_d_cspl(self%nr, self%eqpsi, self%pressure, 1, [rp], p_eqpsi, dp_deqpsi)
       do i = -nth_used, nth_used
          grad(i, 1:2) = grad(i, 1:2) * dp_deqpsi(1) * 0.5 * self%beta_0
       end do
    end if
  end subroutine gradient