betapol.f90

Go to the documentation of this file.

subroutine betapol(a, bpl)
!----------------------------------------------------------------------
! subroutine to evaluate output quantities
!----------------------------------------------------------------------
  
  use itm_types
  use mod_parameter
  use mod_corners
  use mod_dat
  use mod_gauss
  use mod_mesh
  use mod_interpolation
  use mod_profiles
  use mod_spline
  use matrix_calculations
  use phys_functions
  use helena_spline

  implicit none
  
  real(r8), intent(in) :: a
  
  real(r8), intent(out) :: bpl
  
  real(r8), dimension(3) :: ablt
  real(r8) :: factas
  real(r8) :: parea, carea
  real(r8) :: r, wr, s, ws
  real(r8) :: x, xr, xs, y, yr, ys, ps, psr, pss
  real(r8) :: xjac
  real(r8) :: arhs
  integer(itm_i4) :: i, j, ngr, ngs
  integer(itm_i4) :: nelm
  integer(itm_i4) :: n1, n2, n3, n4
  
  factas = 2._r8
  if (ias == 1) factas = 1._r8
  
  parea = 0._r8
  carea = 0._r8

!-- spline p_int and gam_int with use of known derivatives
  call allocate_spline_coefficients(p_spline, npts)
  call spline(npts, psivec, p_int, dpres(1) * sign(1._r8, cpsurfin), &
   dpres(npts) * sign(1._r8, cpsurfin), 1, p_spline)

!------------------------------------- nelm elements ------------------
  do i = nr - 1, 1, -1
    do j = 1, np - 1
      nelm = (i - 1) * (np - 1) + j
      call set_node_number(nelm, n1, n2, n3, n4)
!------------------------------------- 4 point gaussian int. in r -----
      do ngr = 1, 4
        r = xgauss(ngr)
        wr = wgauss(ngr)
!------------------------------------- 4 point gaussian int. in s -----
        do ngs = 1, 4
          s  = xgauss(ngs)
          ws = wgauss(ngs)
          call interpolation(1, xx(1, n1), xx(1, n2), xx(1, n3), xx(1, n4), &
           r, s, x, xr, xs)
          call interpolation(1, yy(1, n1), yy(1, n2), yy(1, n3), yy(1, n4), &
           r, s, y, yr, ys)
          call interpolation(1, 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)
          xjac = xr * ys - xs * yr     
          parea = parea + wr * ws * xjac * spwert(npts, ps, p_spline, &
           psivec, ablt, 0)
          arhs = rh_side(a, x, 0._r8, 0._r8, 0._r8, 0._r8, ps, 0._r8, 0._r8)
          carea = carea + wr * ws * xjac * arhs
        end do
      end do
    end do
  end do
  if (hbt) then
    parea = 0.5_r8 * factas * a * b * parea
    carea = factas * a * carea
  else
    parea = factas * eps * a * b * parea
    carea = factas * a * carea
  end if
  bpl = -(8._r8 * pi / eps) * parea / carea**2

  call deallocate_spline_coefficients(p_spline)

  return
end subroutine betapol