FIXME : Add documentation | Designed for eeq/gs2deq only
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real, | intent(in), | dimension(:) | :: | r | ||
real, | intent(in), | dimension(:) | :: | z | ||
real, | intent(in) | :: | Z_mag | |||
real, | intent(out) | :: | a |
subroutine a_minor(r, z, Z_mag, a)
use mp, only: mp_abort
use splines, only: new_spline, splint, delete_spline, spline
implicit none
real, dimension(:), intent (in) :: r, z
real, intent(in) :: Z_mag
real, intent(out) :: a
real :: r1, r2
integer, parameter :: nz = 5, half_nz = int(nz/2.0), invalidPoint = -nz
real, dimension(nz) :: rtmp, ztmp
integer :: i, j, i1, n, k, index
type (spline) :: spl
!CMR, 28/10/08: add code to avoid duplicate points in 5 point spline
! to determine r2 on inboard mid-plane
logical, parameter :: debug=.false.
k = 0 ; n = size(r)
if (debug) write(6,*) "aminor:"
if (debug) write(6,fmt='("r=",10f8.4)') r
if (debug) write(6,fmt='("z=",10f8.4)') z
if (n < nz) call mp_abort('nbbbs < nz -- very strange. Stopping. Look in geo_utils.fpp', .true.)
j = 0
do i = half_nz + 1, 1, -1
j = j + 1
ztmp(j) = z(i) ; rtmp(j) = r(i)
end do
if (r(1) == r(n) .and. z(1) == z(n)) then
do i = n - 1, n - half_nz, -1
j = j + 1
ztmp(j) = z(i) ; rtmp(j) = r(i)
end do
else
do i = n, n - half_nz + 1, -1
j = j + 1
ztmp(j) = z(i) ; rtmp(j) = r(i)
end do
end if
if (debug) write(6,fmt='("rtmp=",5f8.4)') rtmp
if (debug) write(6,fmt='("ztmp=",5f8.4)') ztmp
if (debug) write(6,fmt='("Z_mag=",f8.4)') Z_mag
spl = new_spline(ztmp, rtmp)
r1 = splint(Z_mag, spl)
call delete_spline(spl)
! find point near magnetic axis elevation on low field side
! Initialise i1 to some invalid value so we can check if an appropriate
! point is not found.
i1 = invalidPoint
! Note there appears to be an assumption here about the order
! of z which may not be true
do i = nz, n
if (z(i) - Z_mag > 0.) then
i1 = i - 1
exit
end if
end do
if (i1 == invalidPoint) then
write(*, '("Could not find point near magnetic axis elevation on low field side in aminor of eeq.f90")')
call mp_abort("Could not find point near magnetic axis elevation on low field side in aminor of eeq.f90")
end if
!CMR, 28/10/08: modify this code to avoid duplicate points in 5 point spline
! to determine r2 on inboard mid-plane
index = i1 - half_nz + k
rtmp(1) = r(index)
ztmp(1) = z(index)
do i = 2, nz
index = i1 - half_nz + i - 1 + k
rtmp(i) = r(index)
ztmp(i) = z(index)
! If we find a duplicate point then we increase the
! duplicate point counter in k, used to offset index and
! store the next value instead Note, here we assume we can't
! get three duplicate values. A better approach may be to
! instead use a while loop where we keep going until we've
! found enough points.
if ((rtmp(i)-rtmp(i-1))**2+(ztmp(i)-ztmp(i-1))**2 < 1.0e-7) then
k = k + 1
if (debug) write(6,fmt='("a_minor: duplicate pt, set k=",i2)') k
rtmp(i) = r(index + 1) ; ztmp(i) = z(index + 1)
end if
end do
if (debug) write(6,fmt='("a_minor: rtmp=",5f8.4)') rtmp
if (debug) write(6,fmt='("a_minor: ztmp=",5f8.4)') ztmp
spl = new_spline(ztmp, rtmp)
r2 = splint(Z_mag, spl)
call delete_spline(spl)
a = (r2 - r1) / 2.
if (debug) write(6,*) "a_minor: r1,r2,a=",r1,r2,a
end subroutine a_minor