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')
allocate(darr, source = self%dbbish)
case('P', 'R')
allocate(darr, source = self%dpbish)
case('T')
allocate(darr, source = self%dtbish)
end select
else
select case(char)
case('B')
allocate(darr, source = self%dbcart)
case('P', 'R')
allocate(darr, source = self%dpcart)
case('T')
allocate(darr, source = 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%eqpsi, self%pressure, [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