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 | Intent | Optional | 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 |
pure subroutine derm(self, f, dfm, char)
use constants, only: pi
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