mercier.f90

Go to the documentation of this file.

subroutine mercier(a, dimerc, drmerc, hh, qq, dq, geonc, zjpar)
!----------------------------------------------------------------------
! subroutine for the ideal and resistive mercier criterion
!----------------------------------------------------------------------

  use itm_types
  use mod_parameter
  use mod_dat
  use mod_gauss
  use mod_mesh
  use mod_helena_io, only : out_he
  use mod_interpolation
  use mod_output

  implicit none

  real(r8), intent(in) :: a

  real(r8), dimension(nr), intent(out) :: dimerc, drmerc, hh, qq, dq
  real(r8), dimension(nr), intent(out) :: geonc
  real(r8), dimension(nr), intent(inout) :: zjpar

  real(r8), dimension(nr) :: p0tmp, ftmp, dptmp, dftmp
  real(r8), dimension(nr) :: zps, zjj1, zjj2, zjj3, zjj4, zjj5, zjj6, djj5
  real(r8), dimension(nr) :: bigg
  real(r8) :: factas
  real(r8) :: r0, b0
  real(r8) :: radius, cpsurf
  real(r8) :: rbscale
  real(r8) :: r, s, ws
  real(r8) :: sumq, sumj1, sumj2, sumj3, sumj4, sumj5, sumj6, sumj5r, &
   sumqr, sumb0, sumbop, sumor2
  real(r8) :: x, xr, xs, xrs, xrr, xss
  real(r8) :: y, yr, ys, yrs, yrr, yss
  real(r8) :: ps, psr, pss, psrs, psrr, psss
  real(r8) :: xjac, xjacr, dl, bigr, gradps2, dssdr
  real(r8) :: zjdchi, b02
  real(r8) :: or2av, b02av, gbar
  real(r8) :: shear
  integer(itm_i4) :: i, j, ngs
  integer(itm_i4) :: n1, n2 ,n3 ,n4
  integer(itm_i4) :: ni

  factas = 2._r8
  if (ias == 1) factas = 1._r8

  r0 = 1._r8
  b0 = 1._r8
  radius = eps * r0
  cpsurf = radius**2 * b0 / alfa

  call profiles(p0tmp, ftmp, dptmp, dftmp, a)

  pscale  = b0**2 * eps / alfa**2
  rbscale = b0 * r0
  do i = 1, nr
    p0tmp(i) = p0tmp(i) * pscale
    dptmp(i) = dptmp(i) * pscale
    ftmp(i) = ftmp(i) * rbscale
    dftmp(i) = dftmp(i) * rbscale
  end do

  if (standard_output) then
    write(out_he, 1)
    write(out_he, 2)
    write(out_he, 1) 
    write(out_he, *) '*i,    ss,     d_i,     d_r,     h,',  &
     '       q,   shear,  geonc    *'
    write(out_he, 1)
  end if
!----------------------------------- flux surface integrals
  r = -1._r8
  do i = 1, nr - 1
    sumq = 0._r8
    sumj1 = 0._r8
    sumj2 = 0._r8
    sumj3 = 0._r8
    sumj4 = 0._r8
    sumj5 = 0._r8
    sumj6 = 0._r8
    sumj5r = 0._r8
    sumqr = 0._r8
    sumb0 = 0._r8
    sumbop = 0._r8
    sumor2 = 0._r8
    do j = 1, np - 1
      n1 = (i - 1) * np + j
      n2 = n1 + 1
      n3 = n2 + np
      n4 = n1 + np
!------------------------------------- 4 point gaussian int. in s -----
      do ngs = 1, 4
        s = xgauss(ngs)
        ws = wgauss(ngs)
        call interpolation(2, xx(1, n1), xx(1, n2), xx(1, n3), xx(1, n4), &
         r, s, x, xr, xs, xrs, xrr, xss)
        call interpolation(2, yy(1, n1), yy(1, n2), yy(1, n3), yy(1, n4), &
         r, s, y, yr, ys, yrs, yrr, yss)
        call interpolation(2, psi(4 * (n1 - 1) + 1), psi(4 * (n2 - 1) + 1), &
         psi(4 * (n3 - 1) + 1), psi(4 * (n4 - 1) + 1), r, s, ps, &
         psr, pss, psrs, psrr, psss)

        xjac = xr * ys - xs * yr
        xjacr = xrr * ys + xr * yrs - xrs * yr - xs * yrr

        dl = sqrt(xs**2 + ys**2)
        bigr = 1._r8 + eps * x
        gradps2 = psr**2 * (xs**2 + ys**2) / xjac**2
        dssdr = abs(psr) / (2._r8 * sqrt(ps))          
!------------------------------------------------ normalisations          
        gradps2 = gradps2 * (cpsurf / radius)**2
        xjac = xjac * radius**2
        bigr = bigr
        dl = radius * dl
        psr = cpsurf * psr
        psrr = cpsurf * psrr
        xjacr = xjacr * radius**2
!---------------------------------------------------------------

        zjdchi = bigr * xjac / abs(psr)

        sumj1 = sumj1 - ws * zjdchi / (gradps2 * bigr**2)
        sumj2 = sumj2 - ws * zjdchi / gradps2
        sumj3 = sumj3 - ws * zjdchi * bigr**2 / gradps2
        sumj4 = sumj4 - ws * zjdchi / bigr**2
        sumj5 = sumj5 - ws * zjdchi
        sumj6 = sumj6 - ws * zjdchi * gradps2 / bigr**2 
        sumq = sumq - ws * xjac / (bigr * abs(psr))

        sumqr = sumqr + psrr * xjac / ((psr**2) * bigr) * ws
        sumqr = sumqr - xjacr / (bigr * psr) * ws
        sumqr = sumqr + xjac * eps * xr / ((bigr**2) * psr) * ws

        sumj5r = sumj5r + xjac * bigr * psrr / (psr**2) * ws
        sumj5r = sumj5r - xjacr * bigr / psr * ws
        sumj5r = sumj5r - xjac * eps * xr / psr * ws

!----------------------------------------- d_nc from Hegna's

        b02 = (ftmp(nr - i + 1) / bigr)**2 + gradps2 / bigr**2

        sumb0 = sumb0 - ws * zjdchi * b02
        sumbop = sumbop - ws * zjdchi * b02 / gradps2
        sumor2 = sumor2 - ws * zjdchi / bigr**2

      end do
    end do
    ni = nr - i + 1
    zps(ni) = sqrt(ps)
    zjj1(ni) = factas * sumj1 / (2._r8 * pi)
    zjj2(ni) = factas * sumj2 / (2._r8 * pi)
    zjj3(ni) = factas * sumj3 / (2._r8 * pi)
    zjj4(ni) = factas * sumj4 / (2._r8 * pi)
    zjj5(ni) = factas * sumj5 / (2._r8 * pi)
    zjj6(ni) = factas * sumj6 / (2._r8 * pi)

    djj5(ni) = factas * sumj5r / dssdr / (2._r8 * pi)

    qq(ni) = factas * ftmp(ni) * sumq / (2._r8 * pi)
    dq(ni) = dftmp(ni) * sumq + ftmp(ni) * sumqr / dssdr
    dq(ni) = dq(ni) * factas / (2._r8 * pi)

    or2av = factas * sumor2 / (2._r8*pi) / zjj5(ni)
    b02av = factas * sumb0  / (2._r8*pi) / zjj5(ni)
    gbar  = factas * sumbop / (2._r8*pi) 

    geonc(ni) = gbar / b02av 

    zjpar(ni) = (-dptmp(ni) * ftmp(ni) - dftmp(ni) * b02av) &
     / (2._r8 * cpsurf * sqrt(ps))

    bigg(ni) = factas * sumbop / (2._r8 * pi)

    dimerc(ni) = (dptmp(ni) * ftmp(ni) * zjj2(ni) / dq(ni) &
     - 0.5_r8)**2 + dptmp(ni) / dq(ni)**2 * (djj5(ni) - dptmp(ni) &
     * zjj3(ni)) * (ftmp(ni)**2 * zjj1(ni) + zjj4(ni))

    hh(ni) = ftmp(ni) * dptmp(ni) / dq(ni) * (zjj2(ni) & 
      - zjj5(ni) * (zjj4(ni) + ftmp(ni)**2 * zjj1(ni)) & 
      / (zjj6(ni) + ftmp(ni)**2 * zjj4(ni)))

    drmerc(ni) = dimerc(ni) - (hh(ni) - 0.5_r8)**2 

    shear = sqrt(ps) * dq(ni) / qq(ni)

    if (standard_output) &
     write(out_he, 3) i, sqrt(ps), -dimerc(ni), -drmerc(ni), hh(ni), &
     qq(ni), shear, geonc(ni)
  end do

  if (standard_output) write(out_he, 1)

! here previously plotted: mercier criteria dimerc(5),drmerc(5),
!                                           hh(5), qq(5)
  
!----------------------- write data for analysis of Hegna's paper to file
  if (additional_output) then
    write(41, *)  ' nr :'
    write(41, *)  nr - 2
    write(41, *) ' di, dr, h, geonc : '
    do i = 2, nr - 1
      write(41, 42) i, dimerc(i), drmerc(i), hh(i), geonc(i)
    end do
  end if

  1 format(' ', 67('*'))
  2 format(' * ideal and resistive mercier criterion ', 26(' '), '*')
  3 format(i3, f10.5, 2es10.2, 3f8.3, es10.2, 2f8.3, es10.2)
 42 format(i4, 12es12.4)

  return
end subroutine mercier