derm Subroutine

public pure subroutine derm(self, f, dfm, char)

Uses

Calculate the derivative of f w.r.t to the radial and poloidal indexes (i.e. calculate the finite differences without dividing by delta psi and delta theta). - dfm(:,:,1) is the psi derivative - dfm(:,:,2) is the theta derivative - char specifies the periodicity in theta - 'E', 'O' mean continuous at theta = +/- pi - 'T' means a jump of 2 pi at theta = +/- pi

Type Bound

abstract_geo_type

Arguments

Type IntentOptional Attributes Name
class(abstract_geo_type), intent(in) :: self
real, intent(in), dimension(:,:) :: f
real, intent(out), dimension(:,:,:) :: dfm
character(len=1), intent(in) :: char

Contents

Source Code


Source Code

  pure subroutine derm(self, f, dfm, char)
    use constants, only: pi, twopi
    implicit none
    class(abstract_geo_type), intent(in) :: self
    real, dimension(:,:), intent(in) :: f
    real, dimension(:,:,:), intent(out) :: dfm
    character(1), intent(in) :: char
    integer :: i, j, nr, nt
    nr = self%nr ; nt = self%nt

    i = 1
    dfm(i,:,1) = -0.5 * (3 * f(i,:) - 4 * f(i+1,:) + f(i+2,:))

    i = nr
    dfm(i,:,1) = 0.5 * (3 * f(i,:) - 4 * f(i-1,:) + f(i-2,:))

    do i = 2, nr-1
       dfm(i,:,1) = 0.5 * (f(i+1,:) - f(i-1,:))
    end do

    do j = 2, nt-1
       dfm(:,j,2) = 0.5 * (f(:,j+1) - f(:,j-1))
    end do

    ! Handle periodic boundary
    if (self%has_full_theta_range) then
       ! Compute df/dtheta at -pi,+pi assuming
       !   (i) periodicity in theta,
       !  (ii) endpoints at -pi,+pi
       ! (iii) equispaced grid in theta, and
       !  (iv) 1 and nt correspond to -pi, +pi and are at the same place
       dfm(:,1,2) = 0.5 * (f(:,2) - f(:,nt-1))
       dfm(:,nt,2) = dfm(:, 1, 2)
       if (char == 'T') then
          dfm(:,1,2) = dfm(:,1,2) + pi
          dfm(:,nt,2) = dfm(:,nt,2) + pi
       end if
    else
       ! assume up-down symmetry for now:
       select case (char)
       case ('E')
          dfm(:,1,2) = 0.0
          dfm(:,nt,2) = 0.0
       case ('O')
          dfm(:,1,2) = f(:,2)
          dfm(:,nt,2) = -f(:,nt-1)
       case ('T')
          dfm(:,1,2) = f(:,2) !f(:,j+1) - f(:,j) -- assumes f(:,1) = 0
          dfm(:,nt,2) = pi - f(:,nt-1) !f(:, j) - f(:, j-1)  -- assumes f(:,nt) = pi
       end select
    end if
  end subroutine derm