suydam_ballooning.f90

Go to the documentation of this file.

subroutine suydam_ballooning(equilibrium_out, zvol, zvolp, xaxis)
!-----------------------------------------------------------------------
! program to determine the ballooning stability of helena equilibria
!         - reads the mapping file as used by castor
!         - calculates stability index using suydam method
!         - for symmetric and asymmetric plasma boundary shapes
!
!         version : 1                 date : 28-09-95
!-----------------------------------------------------------------------
! status : 
!
! 28/9/95   - tested for symmetric soloviev equilibrium with e=1.41
!             eps=0.382, compares well with turnbull results.
!           - also tested for sym. soloviev asymmetric helena
!
!-----------------------------------------------------------------------

  use itm_types
  use mod_parameter
  use mod_dat
  use mod_map
  use mod_mesh
  use mod_profiles
  use mod_spline
  use mod_suydam
  use mod_helena_io, only : out_he
  use mod_output

  use euitm_schemas

  implicit none

  type(type_equilibrium), intent(in) :: equilibrium_out

  real(r8), dimension(nr), intent(in) :: zvol, zvolp
  real(r8), intent(in) :: xaxis

  integer(itm_i4) :: i, j
  real(r8) :: delf, fact
  real(r8) :: funst, fstab, fmarg
  real(r8) :: shear, alpha, shear1, alpha2
!  real :: alpha1, alpha11
  real(r8) :: rho, drho_ds
  real(r8) :: tbb, tbf
  character(len = 25) :: bal,inibal

  call suydam_initialization(ias, tbb, tbf)

!---------------------- read helena mapping file -----------------------      
  call suydam_read_helena(equilibrium_out)

!---------------------- calculate b-field on grid points ---------------     
  call suydam_bfield

!---------------------- calculate p and q coeff. on all gridpoints -----      
  call suydam_pq(ias)

!---------------------- calculate stability index ----------------------     
  if (standard_output) then
    write(out_he,*) '****************************************************', &
     '************************************'
    write(out_he,*) '* i,   flux,   rho,    q,    shear,   shear1,', &
     '   alpha,   alpha1,   fmarg,    ballooning *'
    write(out_he,*) '****************************************************', &
     '************************************'
  end if
  do i = 2, npsi
    fact = 1._r8
    call suydam(i, 0._r8, tbb, tbf, 1._r8, bal)
!---------------------- find distance from stability boundary ----------   
    inibal = bal
    if (bal == ' stable') then
      delf = 2._r8
    else
      delf = 0.5_r8
    end if
    fact = delf * fact
    do j = 1, 10
      call suydam(i, 0._r8, tbb, tbf, fact, bal)
      if (bal == inibal) then
        fact = delf * fact
      else
        goto 10
      end if
    end do
!--------------------- upper and lower limits established --------------    
10  if (inibal == ' stable') then
      funst = fact
      fstab = fact / 2._r8
    else
      fstab = fact
      funst = fact * 2._r8
    end if   
!------------------------------------- bisection to find factor
    do j = 1, 10 
      fact = (funst + fstab) / 2._r8
      call suydam(i, 0._r8, tbb, tbf, fact, bal)
      if (bal == ' stable') then
        fstab = fact
      else
        funst = fact
      end if               
    end do      
    shear = cs(i) / qs(i) * dqs(i) 
    alpha = -2._r8 * qs(i)**2 * p_spline%sp2(i) / eps

    shear1 = 2._r8* zvol(i) / qs(i) * dqs(i) / (2._r8 * cs(i) *zvolp(i))
!-------------------------------------------- lao's definition
!       alpha1 = - p_spline%sp2(i) / (2._r8 * cs(i)) * zvolp(i) / cpsurf**2
!     >  * sqrt(zvol(i) / (4._r8 * pi * (1._r8 + eps *xaxis))) * eps**3.
!
!--------------------------------------------- use rho as radius
    rho = sqrt(zvol(i) / zvol(nr))
    drho_ds =  cs(i) / (rho * zvol(nr)) * zvolp(i)  
    alpha2 = alpha / drho_ds
!----------------------------------------- lao corrected? needs check
!       alpha11 = - p_spline%sp2(i) / (2_r8 * cs(i)) * zvolp(i) / cpsurf**2
!     >  / (pi**2 * eps) * rho * zvol(nr) * eps**4
!----------------------------------------------------------------------

    fmarg = (fstab + funst) / 2._r8

!-----------------------------temporary fix for negative shear and ias=1

    if (shear < 0._r8 .and. ias == 1) then
      fmarg = 0._r8
      inibal =  ' stable'
    end if

    if (standard_output) &
     write(out_he,11) i, cs(i)**2, rho, qs(i), shear, shear1, &
     alpha, alpha2, fmarg, inibal
  end do

  if (standard_output) &
   write(out_he, *) '****************************************************', &
   '*****************'

  deallocate(b02)
  deallocate(dt_b02)
  deallocate(ds_b02)
  deallocate(cp0)
  deallocate(cp1)
  deallocate(cp2)
  deallocate(cq0)
  deallocate(cq1)

  call deallocate_spline_coefficients(q_spline)
  call deallocate_spline_coefficients(p_spline)
  call deallocate_spline_coefficients(rbphi_spline)

 11 format(' ',i3,' ',f9.6,' ',f9.6,' ',f7.4,' ',f8.4,' ',f8.4,' ', &
         e10.4,' ',e10.4,' ',f8.3,'  ',a25)

  return
end subroutine suydam_ballooning