check_stability Subroutine

  subroutine check_stability(shat_in, beta_prime_in, dbetadrho_out, &
       iunstable, theta0, restore)
    !Note iunstable>0 means unstable - Could use a logical but size may be interesting
    use theta_grid_params, only: ntheta, nperiod
    use geometry, only: eikcoefs, eikcoefs_output_type
    use optionals, only: get_option_with_default
    implicit none
    real, intent(in) :: shat_in, beta_prime_in
    real, intent(out) :: dbetadrho_out
    integer, intent(out) :: iunstable
    real, intent(in), optional :: theta0
    logical, intent(in), optional :: restore
    type(eikcoefs_output_type) :: eikcoefs_results
    real :: shat_bak, beta_prime_bak
    integer :: ntgrid
    logical :: local_restore
    
    local_restore = get_option_with_default(restore, .false.)

    !Store initial state if we're going to restore
    if (local_restore) then
       call get_shat(shat_bak)
       call get_beta_prime(beta_prime_bak)
    end if

    !Setup geometry
    !NOTE: As how we can influence beta_prime (and shat) varies with bishop
    !we actually end up setting different parameters in different cases. The
    !physical parameter dbetadrho should really be the relevant parameter.
    call set_beta_prime(beta_prime_in)
    call set_shat(shat_in)
    call eikcoefs(ntheta, nperiod, eikcoefs_results) !Recalculate geometry terms

    !Here we pass out the value of dbetadrho, this should be bishop value independent
    dbetadrho_out = eikcoefs_results%dbetadrho

    !Note we don't use theta_grid:ntgrid as in the case that
    !gridgen resizes the theta grid used by GS2 we don't end up
    !resizing any of the geometry module arrays (that we use directly here)
    !Really we would like to use everything from theta_grid directly but we'd
    !need some way to force theta_grid to update it's internal arrays when we
    !change parameters
    ntgrid=ubound(eikcoefs_results%theta,1)

    !Test stability
    iunstable = is_unstable(ntgrid, eikcoefs_results%theta, eikcoefs_results%bmag, &
         eikcoefs_results%dbetadrho, eikcoefs_results%gradpar, eikcoefs_results%gds2, &
         eikcoefs_results%gds21, eikcoefs_results%gds22, eikcoefs_results%cvdrift, &
         eikcoefs_results%cvdrift0, theta0)

    !Restore geometry if requested
    if (local_restore) then
       call set_beta_prime(beta_prime_bak)
       call set_shat(shat_bak)
       call eikcoefs(ntheta, nperiod, eikcoefs_results)
    end if
  end subroutine check_stability