flxint2.f90

Go to the documentation of this file.

subroutine flxint2(xaxis, a, cx, cy, sumdq)
!----------------------------------------------------------------------
! subroutine to evaluate new fdf profile
!----------------------------------------------------------------------

  use itm_types
  use mod_parameter
  use mod_dat
  use mod_gauss
  use mod_mesh
  use mod_profiles
  use mod_interpolation
  use helena_spline

  implicit none

  real(r8), intent(in) :: xaxis, a, cx, cy

  real(r8), intent(out) :: sumdq

  type(spline_coefficients) :: dd_spline

  real(r8), dimension(nr) :: qplin, deltaq, df2old
  real(r8), dimension(nr) :: qplot, xplot, df2tmp
  real(r8), dimension(nr) :: zps
  real(r8), dimension(3) :: abltg
  real(r8) :: factas
  real(r8) :: r, s, ws
  real(r8) :: sumq, sumqr
  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, er
  real(r8) :: qtmp, qtmp0
  real(r8) :: ss, dummy, df21
  integer(itm_i4) :: i, j, ngs
  integer(itm_i4) :: n1, n2, n3, n4

  real(r8) :: q_profile, dgamma_dpsi

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

  call allocate_spline_coefficients(dd_spline, nr)
  call profiles(dd_spline%sp1, dd_spline%sp2, dd_spline%sp3, &
   dd_spline%sp4, a)
!------------------------------------- nelm elements ------------------
  r = -1._r8
  do i = 1, nr - 1
    sumq = 0._r8
    sumqr = 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
        bigr = 1._r8 + eps * x
        sumq = sumq - ws * xjac / ( bigr * abs(psr))
        er = xrr * ys + xr * yrs - xrs * yr - xs * yrr
        sumqr = sumqr + psrr * xjac / ((psr**2) * bigr) * ws
        sumqr = sumqr - er / (bigr * psr) * ws
        sumqr = sumqr + xjac * eps * xr / ((bigr**2) * psr) * ws
      end do
    end do
    zps(nr - i + 1) = ps
    qtmp = 0.5_r8 * factas * (sumq * alfa / pi) * dd_spline%sp2(nr - i + 1)
    qplin(nr - i + 1) = q_profile(ps)
    deltaq(nr - i + 1) = qplin(nr - i + 1) - qtmp 
    df2old(nr - i + 1) = dgamma_dpsi(ps)
    qplot(nr - i + 1) = qtmp
    xplot(nr - i + 1) = sqrt(ps)
  end do
  qtmp0 = dd_spline%sp2(1) * (alfa / (2._r8 * sqrt(cx * cy) * (1._r8 &
   + eps * xaxis)))
  qplot(1) = qtmp0
  xplot(1) = 0._r8
  zps(1) = 0._r8
  qplin(1) = q_profile(0._r8)
  deltaq(1) = qplin(1) - qtmp0
  df2old(1) = 1._r8

  sumdq = 0._r8
  df2tmp(1) = df2old(1) + 0.25_r8 * ampl * (deltaq(1) - deltaq(2))
  zps(1) = 0._r8 
  df2tmp(nr) = df2old(nr)
  do i = 2, nr - 1
    df2tmp(nr - i + 1) = df2old(nr - i + 1) + ampl * (deltaq(nr - i) &
     - deltaq(nr - i + 2)) + ampl * deltaq(nr - i + 1) / 4._r8
    df2tmp(nr - i + 1) = df2tmp(nr - i + 1) / df2tmp(1)
!       write(20, 41) xplot(nr - i + 1), qplot(nr - i + 1), deltaq(nr
!     >   - i + 1), deltaq(nr - i) - deltaq(nr - i + 2), df2tmp(nr - i
!     >   + 1), df2old(nr - i + 1)
    sumdq = sumdq + abs(deltaq(nr - i + 1))
  end do
  df2tmp(1) = 1._r8
  df2tmp(nr - 1) = df2old(nr - 1) + 0.25_r8 * ampl * (deltaq(nr - 2) &
   - deltaq(nr - 1))

  sumdq = sumdq / dble(nr - 2)

!      write(20, 41) zps(1), qtmp0, q_profile(0.), df2tmp(1),
!     >  dgamma_dpsi(0.)

!-------------------------------------- still freedom to scale q-profile
!      alfa = alfa * qplin(1) / qplot(1)
!      write(20, *) alfa, qplin(1), qplot(1)
!------------------------------------- calc equidistant gamma profile   
  call spline(nr, xplot, df2tmp, 0._r8, 0._r8, 2, dd_spline)
  do i = 1, npts
    ps = dble(i - 1) / dble(npts - 1)
    ss = sqrt(ps)
    dummy = spwert(nr, ss, dd_spline, xplot, abltg, 0)
    fdf(i) = dummy
!    write(20, 41) ps, ss, fdf(i), df2tmp(i)
  end do
  df21 = fdf(1)
  do i = 1, npts
    fdf(i) = fdf(i) / df21
  end do

  call deallocate_spline_coefficients(dd_spline)

! 41 format(12e12.4)

  return
end subroutine flxint2