mhdballn Subroutine

  subroutine mhdballn(ntgrid, beta_prime, diffscheme, iunstable, psiend, cflmax, &
       cflmin, psi, bmag, cvdrift, gradpar, gds2, theta)
    implicit none
    integer, intent(in) :: ntgrid
    real, intent(in) :: beta_prime, diffscheme
    integer, intent(out) :: iunstable
    real, intent(out) :: psiend, cflmax, cflmin
    real, dimension(-ntgrid:ntgrid), intent(out) :: psi
    real, dimension(-ntgrid:ntgrid), intent(in) :: bmag, cvdrift, gradpar, gds2, theta
    integer :: i
    real, dimension(-ntgrid:ntgrid) :: psi_p, c, g, c1, c2, ch, g1, g2, gh, delx
    real :: alpha,  b_prime

    !       include 'ball.inc'

    !       input arguments:
    !       npt=total number of gridpoints with theta>0 = nperiod*ntheta+1
    !       beta_prime=d beta/ d rho
    !       alpha=numerical knob (0 to .333 are good if cflmax,min<1)
    !
    !       output arguments:
    !       iunstable=number of times the solution crossed the axis. If
    !                 this is greater than 0 it is unstable.
    !
    !       psi is the solution of the Euler Lagrange equation for delta W.
    !       At the first theta point, psi is zero with slope positive 1.
    !
    !       psiend=psi at the last gridpoint
    !       psi is psi at every grid point
    !       cflmax,min are indicators of whether resolution is adequate.
    !       diffscheme is a numerical difference scheme knob.

    !       For greatest accuracy diffscheme=0. or .333333 .
    !       If the absolute value of cflmax,min approaches one, inaccuracy and/or
    !       numerical instability are possible.  If cflmax,min are >1-3, increase
    !       diffscheme to .666 or 1.0 for stability, but probably resolution is
    !       inadequate.
    !

    b_prime = -beta_prime
    iunstable=0
    alpha=1.-diffscheme

    g = gds2*gradpar/bmag
    c = b_prime*cvdrift/(2.*bmag*gradpar)

    i=-ntgrid
    cflmax=(theta(i+1)-theta(i))**2*(c(i+1)+c(i))/(g(i+1)+g(i))
    cflmin=cflmax
    do i=-ntgrid+1, ntgrid
       delx(i)=theta(i)-theta(i-1)
       ch(i)=.5*(c(i)+c(i-1))
       gh(i)=.5*(g(i)+g(i-1))
       cflmax=max(cflmax,delx(i)**2*ch(i)/gh(i))
       cflmin=min(cflmin,delx(i)**2*ch(i)/gh(i))
    enddo


    do i=-ntgrid+1, ntgrid-1
       c1(i)=-delx(i+1)*(alpha*c(i)+.5*diffscheme*ch(i+1)) &
            -delx(i  )*(alpha*c(i)+.5*diffscheme*ch(i  )) &
            -delx(i)*.5*diffscheme*ch(i)
       c1(i)=.5*c1(i)
    enddo

    do i=-ntgrid+1, ntgrid
       c2(i)=-.25*diffscheme*ch(i)*delx(i)
       g1(i)=gh(i)/delx(i)
       g2(i)=1./(gh(i)/delx(i)+.25*diffscheme*ch(i)*delx(i))
    enddo

    i=-ntgrid
    psi(i)=0.
    i=i+1
    psi(i)=delx(i)
    psi_p(i)=psi(i)/g2(i)-g1(i)*psi(i-1)

    do i=-ntgrid+1,ntgrid-1
       psi_p(i+1)=psi_p(i)+c1(i)*psi(i)+c2(i)*psi(i-1)
       psi(i+1)=(g1(i+1)*psi(i)+psi_p(i+1))*g2(i+1)
    enddo
    psiend=psi(ntgrid)

    do i=-ntgrid+1,ntgrid-1
       if(psi(i)*psi(i+1) <= 0.) iunstable=iunstable+1
    enddo

  end subroutine mhdballn