is_unstable Function

  pure integer function is_unstable(ntgrid, theta, bmag, dbetadrho, gradpar, &
       gds2, gds21, gds22, cvdrift, cvdrift0, theta0_in) &
       result(is_unstable_out)
    use optionals, only: get_option_with_default
    implicit none
    integer, intent(in) :: ntgrid
    real, dimension(-ntgrid:ntgrid), intent(in) :: theta, bmag, gradpar
    real, intent(in) :: dbetadrho
    real, dimension(-ntgrid:ntgrid), intent(in) :: gds2, gds21, gds22, cvdrift, cvdrift0
    real, intent(in), optional :: theta0_in
    !<DD> NEED TO IMPROVE VARIABLE NAMES
    real, dimension(-ntgrid:ntgrid) :: g, c, psi, c1, c2, ch, g1, g2, gh
    real, dimension(-ntgrid:ntgrid-1) :: delthet
    integer :: ig
    real :: one_m_diff, psi_prime, theta0

    !Initialise
    theta0 = get_option_with_default(theta0_in, 0.0)
    is_unstable_out = 0
    one_m_diff = 1 - diff

    !Calculate geometry arrays
    delthet = theta(-ntgrid+1:) - theta(:ntgrid-1)

    g = (gds2 + 2 * theta0 * gds21 + theta0**2 * gds22) * gradpar / bmag
    c = -0.5 * dbetadrho * (cvdrift + theta0 * cvdrift0) / (bmag * gradpar)

    ! Get grid centre values -- note we don't set the value at -ntgrid
    do ig=-ntgrid+1, ntgrid
       ch(ig)=0.5*(c(ig)+c(ig-1))
       gh(ig)=0.5*(g(ig)+g(ig-1))
    enddo

    !Calculate coefficients -- note we don't set the value at -ntgrid or ntgrid
    do ig=-ntgrid+1,ntgrid-1
       c1(ig) = -delthet(ig)*(one_m_diff*c(ig)+0.5*diff*ch(ig+1)) &
            -delthet(ig-1)*(one_m_diff*c(ig)+0.5*diff*ch(ig)) &
            -delthet(ig-1)*0.5*diff*ch(ig)
       c1(ig)=0.5*c1(ig)
    enddo

    !Calculate coefficients -- note we don't set the value at -ntgrid
    do ig=-ntgrid+1,ntgrid
       c2(ig) = -0.25*diff*ch(ig)*delthet(ig-1)
       g1(ig) = gh(ig)/delthet(ig-1)
       g2(ig) = 1.0/(0.25*diff*ch(ig)*delthet(ig-1)+gh(ig)/delthet(ig-1))
    enddo

    !Set boundary values
    psi(-ntgrid)=0.
    psi(-ntgrid+1)=delthet(-ntgrid)
    psi_prime=psi(-ntgrid+1)/g2(-ntgrid+1)!-g1(-ntgrid+1)*psi(-ntgrid) !Second term zero by defn

    !Integrate along in theta
    do ig=-ntgrid+1,ntgrid-1
       psi_prime=psi_prime+c1(ig)*psi(ig)+c2(ig)*psi(ig-1)
       psi(ig+1)=(g1(ig+1)*psi(ig)+psi_prime)*g2(ig+1)
    enddo

    !Find how many times psi crosses axis || Actually how many times it hits zero
    !If we really do just want to count crossings then need to replace le with lt
    ! Could we replace this with
    ! is_unstable_out = count(psi(-ntgrid+1:ntgrid-1) * psi(-ntgrid+2:) <= 0)
    do ig=-ntgrid+1,ntgrid-1
       if(psi(ig)*psi(ig+1) .le. 0 ) is_unstable_out = is_unstable_out + 1
    enddo

    !Probably want to think about writing out psi if debug
  end function is_unstable