eqitem_eqfile Function

public function eqitem_eqfile(self, r_in, theta_in, f, char) result(fstar)

Uses

Calculates fstar which is f interpolated at the location (r,theta). Here r is the normalised poloidal flux coordinate rp (= psi_pN + psi_0N) and theta_in is theta. f is a grid of values of f as a function of psi_p,theta First we find the radial and theta indices bounding our requested point and we then interpolate within the bound rectangle.

Type Bound

abstract_eqfile_geo_type

Arguments

Type IntentOptional Attributes Name
class(abstract_eqfile_geo_type), intent(in) :: self
real, intent(in) :: r_in
real, intent(in) :: theta_in
real, intent(in), dimension(:,:) :: f
character(len=1), intent(in) :: char

Return Value real


Contents

Source Code


Source Code

  real function eqitem_eqfile(self, r_in, theta_in, f, char) result(fstar)
    use mp, only: mp_abort
    implicit none
    class(abstract_eqfile_geo_type), intent(in) :: self
    real, intent(in) :: r_in, theta_in
    real, dimension(:,:), intent(in) :: f
    character(1), intent(in) :: char
    integer :: istar, jstar
    real :: r, thet, eqpsi_sign, tps
    ! allow psi(r) to be a decreasing function
    eqpsi_sign = 1.
    if (self%eqpsi(2) < self%eqpsi(1)) eqpsi_sign = -1.
    associate (dim_1 => eqpsi_sign * self%eqpsi , mtheta => self%eqth(1, :)) !Non-contiguous copy - slow?

      !Copy input r to avoid any modification
      r = r_in ; thet = mod2pi(theta_in) ; tps = 1.

      if (.not. self%has_full_theta_range) then ! assume up-down symmetry
         if (char == 'Z' .and. abs(thet) >= 1.e-10) tps = sign(1.0, thet)
         thet = abs(thet)
      end if

      ! check for axis evaluation, no gradients available at axis, so do not get too close
      if (r < dim_1(2)) then
         write(*,*) r, theta_in, eqpsi_sign, dim_1(1), dim_1(2)
         call mp_abort('no evaluation too close to or inside axis allowed in eqitem', .true.)
      end if

      ! disallow evaluations on or outside the surface for now
      if (r >= dim_1(self%nr)) then
         write(*,*) r, theta_in, dim_1(self%nr), eqpsi_sign
         call mp_abort('No evaluation of eqitem allowed on or outside surface', .true.)
      end if

      istar = find_lower_bound(r, dim_1, 0)
      jstar = find_lower_bound(thet, mtheta, self%nt - 1)
      fstar = tps * bilinear_interpolate(f, dim_1, mtheta, istar, jstar, r, thet)
    end associate
  end function eqitem_eqfile