passing_particles.f90

Go to the documentation of this file.

subroutine passing_particles(a, fcirc, b02av, b0max, rav, equilibrium)
!----------------------------------------------------------------------
! subroutine for the fraction of circulating particles
!----------------------------------------------------------------------

  use itm_types
  use mod_parameter
  use mod_dat
  use mod_gauss
  use mod_mesh
  use mod_helena_io, only : iu6
  use mod_interpolation
  use euitm_schemas

  implicit none

  type(type_equilibrium), intent(inout) :: equilibrium
  real(r8), intent(in) :: a
  real(r8), dimension(nr), intent(out) :: fcirc, b02av, b0max, rav

  integer(itm_i4), parameter :: nk = 101

  real(r8), dimension(nk) :: sumk
  real(r8), dimension(nr) :: p0tmp, ftmp, dptmp, dftmp
  real(r8) :: factas
  real(r8) :: r0, b0
  real(r8) :: radius, cpsurf
  real(r8) :: rbscale
  real(r8) :: r, s, ws
  real(r8) :: sum1, sum2, sumr, sum
  real(r8) :: b02max, b02, bm
  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, bigr, gradps2, zjdchi, zlam
  integer(itm_i4) :: i, j, k, 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)

  rbscale = b0 * r0

  do i = 1, nr
    ftmp(i) = ftmp(i) * rbscale
    dftmp(i) = dftmp(i) * rbscale
  end do
!---------------------- find bmax and b^2 average on every surface
  r = -1._r8
  do i = 1, nr - 1
    sum1 = 0._r8
    sum2 = 0._r8
    sumr = 0._r8
    b02max = -1.e20_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
        bigr = 1._r8 + eps * x
        gradps2 = psr**2 * (xs**2 + ys**2) / xjac**2
!------------------------------------------------ normalisations          
        gradps2 = gradps2 * (cpsurf / radius)**2
        xjac = xjac * radius**2
        bigr = bigr
        psr = cpsurf * psr
!---------------------------------------------------------------
        zjdchi = bigr * xjac / abs(psr)
        b02 = (ftmp(nr - i + 1) / bigr)**2 + gradps2 / bigr**2
        sum1  = sum1 - ws * zjdchi
        sum2  = sum2 - ws * zjdchi * b02
        sumr  = sumr - ws * zjdchi * bigr

        if (b02 > b02max) b02max = b02
      end do
    end do
    ni = nr - i + 1
    sum1 = factas * sum1 / (2._r8 * pi)
    sum2 = factas * sum2 / (2._r8 * pi)  
    sumr = factas * sumr / (2._r8 * pi)
    b02av(ni) = sum2 / sum1
    b0max(ni) = sqrt(abs(b02max))
    rav(ni) = sumr / sum1
  end do

!---------------------------------- calculate average term in integral
  r = -1._r8
  do i = 1, nr - 1
    sum1 = 0._r8
    do k = 1, nk
      sumk(k) = 0._r8
    end do
    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
        bigr = 1._r8 + eps * x
        gradps2 = psr**2 * (xs**2 + ys**2) / xjac**2
!------------------------------------------------ normalisations          
        gradps2 = gradps2 * (cpsurf / radius)**2
        xjac = xjac * radius**2
        bigr = bigr
        psr = cpsurf * psr
!---------------------------------------------------------------
        zjdchi = bigr * xjac / abs(psr)
        b02 = (ftmp(nr - i + 1) / bigr)**2 + gradps2 / bigr**2
        b0 = sqrt(b02)
        bm = b0max(nr - i + 1)
        do k = 1, nk 
          zlam = dble(k - 1) / dble(nk - 1)
          sumk(k) = sumk(k) - ws * zjdchi * sqrt(abs(1._r8 - zlam &
           * b0 / bm))
        end do
      sum1  = sum1 - ws * zjdchi
      end do
    end do
!--- axis: ni = 1, boundary: ni = nr
    ni = nr - i + 1

!------------------------------------------ integrate over lambda
    sum = 0._r8
    do k = 2, nk - 1
      zlam = dble(k - 1) / dble(nk - 1)
      sum = sum + zlam * sum1 / sumk(k)
    end do
    fcirc(ni) = (sum + 0.5_r8 * sum1 / sumk(nk)) / dble(nk - 1)
    fcirc(ni) = 0.75_r8 * fcirc(ni) * b02av(ni) / b0max(ni)**2
!   if (verbosity > 3) &
!    write(iu6, 41) i, sqrt(ps), b02av(ni), b0max(ni), fcirc(ni)
  end do

  equilibrium%profiles_1d%ftrap(2 : nr) = 1.0_r8 - fcirc(2 : nr)
!--- no trapped particles on axis
  equilibrium%profiles_1d%ftrap(1) = 0._r8

  if (verbosity > 2) write(iu6, *) ' done passing particles'

! 41 format(i3, 15es11.3)

  return
end subroutine passing_particles