diff --git a/Sources/check3d.f90 b/Sources/check3d.f90
index 404b59c..90095ce 100644
--- a/Sources/check3d.f90
+++ b/Sources/check3d.f90
@@ -323,27 +323,27 @@ SUBROUTINE CheckForces(xc, gc)
DO ntype = 1, ndims
DO n = -ntor, ntor
DO m = 0, mpol
- icol = icol+1
-
- xc(icol,js) = eps
- CALL funct_island
- wplus = wtotal
-
- xc(icol,js) = -eps
- CALL funct_island
- wmins = wtotal
- force = -(wplus-wmins)/(2*eps)
-
- CALL ASSERT(l_getwmhd,'L_GETWMHD = FALSE!')
-
- xc(icol,js) = 0
-
- IF (lWrite) THEN
- nt1 = ntype
- IF (ntype .GT. 3) nt1 = ntype-3
- WRITE (3000+iam,100) m, n, ntype, force, force_array(nt1), &
- gc1(icol,js)
- END IF
+ icol = icol+1
+
+ xc(icol,js) = eps
+ CALL funct_island
+ wplus = wtotal
+
+ xc(icol,js) = -eps
+ CALL funct_island
+ wmins = wtotal
+ force = -(wplus-wmins)/(2*eps)
+
+ CALL ASSERT(l_getwmhd,'L_GETWMHD = FALSE!')
+
+ xc(icol,js) = 0
+
+ IF (lWrite) THEN
+ nt1 = ntype
+ IF (ntype .GT. 3) nt1 = ntype-3
+ WRITE (3000+iam,100) m, n, ntype, force, force_array(nt1), &
+ gc1(icol,js)
+ END IF
END DO
END DO
diff --git a/Sources/diagnostics_mod.f90 b/Sources/diagnostics_mod.f90
index 7c05987..c84824a 100644
--- a/Sources/diagnostics_mod.f90
+++ b/Sources/diagnostics_mod.f90
@@ -94,25 +94,6 @@ SUBROUTINE divb(ns_min, ns_max)
f_cos, nsmin, nsmax)
END IF
-! DO m = 0, mpol
-! DO n = -ntor, ntor
-! divbmnsf(m,n,nl:nh) = ohs*(jbsupsmnsh(m,n,nl+1:nh+1) - &
-! jbsupsmnsh(m,n,nl:nh)) &
-! - m*(jbsupumnch(m,n,nl+1:nh+1) + &
-! jbsupumnch(m,n,nl:nh))*0.5 &
-! - n*nfp*(jbsupvmnch(m,n,nl+1:nh+1) + &
-! jbsupvmnch(m,n,nl:nh))*0.5
-! IF (lasym) THEN
-! divbmncf(m,n,nl:nh) = ohs*(jbsupsmnch(m,n,nl+1:nh+1) - &
-! jbsupsmnch(m,n,nl:nh)) &
-! - m*(jbsupumnsh(m,n,nl+1:nh+1) + &
-! jbsupumnsh(m,n,nl:nh))*0.5 &
-! + n*nfp*(jbsupvmnsh(m,n,nl+1:nh+1) + &
-! jbsupvmnsh(m,n,nl:nh))*0.5
-! END IF
-! END DO
-! END DO
-
tnorm = hs_i*SUM(jbsupumnch(:,:,nsmin:nsmax)**2 &
+ jbsupvmnch(:,:,nsmin:nsmax)**2 &
+ jbsupsmnsh(:,:,nsmin:nsmax)**2)
@@ -144,7 +125,7 @@ SUBROUTINE divb(ns_min, ns_max)
!> @param[in] ns_max Upper radial index.
!-------------------------------------------------------------------------------
FUNCTION divdb(ns_min, ns_max)
- USE island_params, ONLY: ohs=>ohs_i
+ USE island_params, ONLY: ohs=>ohs_i, fourier_context
USE hessian, ONLY: gather_array
! Declare Arguments
@@ -156,6 +137,7 @@ FUNCTION divdb(ns_min, ns_max)
INTEGER :: istat
INTEGER :: m
INTEGER :: n
+ INTEGER :: n_mode
INTEGER :: nl
INTEGER :: nh
REAL (dp) :: tnorm
@@ -184,21 +166,22 @@ FUNCTION divdb(ns_min, ns_max)
divbmnsf = 0
divbmncf = 0
- DO m = 0, mpol
- DO n = -ntor, ntor
- divbmnsf(m,n,nl:nh) = -m*(djbsupumnch(m,n,nl+1:nh+1) + &
- djbsupumnch(m,n,nl:nh))*0.5 &
- - n*nfp*(djbsupvmnch(m,n,nl+1:nh+1) + &
- djbsupvmnch(m,n,nl:nh))*0.5 &
- + ohs*(djbsupsmnsh(m,n,nl+1:nh+1) - &
- djbsupsmnsh(m,n,nl:nh))
+ DO n = -ntor, ntor
+ n_mode = fourier_context%tor_modes(n)*nfp
+ DO m = 0, mpol
+ divbmnsf(m,n,nl:nh) = -m*(djbsupumnch(m,n,nl+1:nh+1) + &
+ & djbsupumnch(m,n,nl:nh))*0.5 &
+ & - n_mode*(djbsupvmnch(m,n,nl+1:nh+1) + &
+ & djbsupvmnch(m,n,nl:nh))*0.5 &
+ & + ohs*(djbsupsmnsh(m,n,nl+1:nh+1) - &
+ & djbsupsmnsh(m,n,nl:nh))
IF (lasym) THEN
- divbmncf(m,n,nl:nh) = m*(djbsupumnsh(m,n,nl+1:nh+1) + &
- djbsupumnsh(m,n,nl:nh))*0.5 &
- + n*nfp*(djbsupvmnsh(m,n,nl+1:nh+1) + &
- djbsupvmnsh(m,n,nl:nh))*0.5 &
- + ohs*(djbsupsmnch(m,n,nl+1:nh+1) - &
- djbsupsmnch(m,n,nl:nh))
+ divbmncf(m,n,nl:nh) = m*(djbsupumnsh(m,n,nl+1:nh+1) + &
+ & djbsupumnsh(m,n,nl:nh))*0.5 &
+ & + n_mode*(djbsupvmnsh(m,n,nl+1:nh+1) + &
+ & djbsupvmnsh(m,n,nl:nh))*0.5 &
+ & + ohs*(djbsupsmnch(m,n,nl+1:nh+1) - &
+ & djbsupsmnch(m,n,nl:nh))
END IF
END DO
END DO
@@ -428,7 +411,6 @@ SUBROUTINE SIESTA_PROFILES (iunit, arr_value, label, fsq_total)
DO n = -ntor, ntor
DO m = 0, mpol
m_array(m,n) = m
- n_array(m,n) = n
END DO
END DO
@@ -437,7 +419,8 @@ SUBROUTINE SIESTA_PROFILES (iunit, arr_value, label, fsq_total)
WRITE (p_format,1000) mnmax
WRITE (iunit,p_format) ' ',' ','MPOL--->', m_array
WRITE (p_format,1000) mnmax
- WRITE (iunit,p_format) 'RADIUS','RMS','NTOR--->', n_array
+ WRITE (iunit,p_format) 'RADIUS','RMS','NTOR--->', &
+ & fourier_context%tor_modes
WRITE (p_format,1002) mnmax
DO js = 2, ns
diff --git a/Sources/fourier.f90 b/Sources/fourier.f90
index 87fa857..335266a 100644
--- a/Sources/fourier.f90
+++ b/Sources/fourier.f90
@@ -1084,8 +1084,7 @@ FUNCTION fourier_get_index(this, n)
DO i = 0, SIGN(UBOUND(this%tor_modes, 1), n), SIGN(1, n)
IF (this%tor_modes(i) .eq. n) THEN
- n = this%tor_modes(i)
-
+ n = i
fourier_get_index = .true.
EXIT
END IF
@@ -1138,7 +1137,6 @@ FUNCTION test_fourier()
four => fourier_class(pol, tor, theta, zeta, nfp, .false., tor_modes)
test_fourier = .true.
- DEALLOCATE(tor_modes)
DO n = -tor, tor
DO m = 0, pol
@@ -1153,7 +1151,7 @@ FUNCTION test_fourier()
CALL four%tomnsp(testij, result, f_cos)
test_fourier = test_fourier .and. &
check(testmn, result, pol, tor, &
- 'cosine_parity_test')
+ & tor_modes, 'cosine_parity_test')
! Test sine parity transform.
IF (m .eq. m0 .and. n .eq. n0) THEN
@@ -1165,7 +1163,7 @@ FUNCTION test_fourier()
CALL four%tomnsp(testij, result, f_sin)
test_fourier = test_fourier .and. &
check(testmn, result, pol, tor, &
- 'sine_parity_test')
+ & tor_modes, 'sine_parity_test')
! Test u derivatives of cosine parity.
testmn = 0
@@ -1175,7 +1173,7 @@ FUNCTION test_fourier()
testmn(m,n) = -m
test_fourier = test_fourier .and. &
check(testmn, result, pol, tor, &
- 'du_cosine_parity_test')
+ & tor_modes, 'du_cosine_parity_test')
! Test v derivatives of cosine parity.
testmn = 0
@@ -1185,7 +1183,7 @@ FUNCTION test_fourier()
testmn(m,n) = -n*nfp
test_fourier = test_fourier .and. &
check(testmn, result, pol, tor, &
- 'dv_cosine_parity_test')
+ & tor_modes, 'dv_cosine_parity_test')
! Test u derivatives of sine parity.
testmn = 0
@@ -1195,7 +1193,7 @@ FUNCTION test_fourier()
testmn(m,n) = m
test_fourier = test_fourier .and. &
check(testmn, result, pol, tor, &
- 'du_sine_parity_test')
+ & tor_modes, 'du_sine_parity_test')
! Test v derivatives of sine parity.
testmn = 0
@@ -1204,8 +1202,8 @@ FUNCTION test_fourier()
CALL four%tomnsp(testij, result, f_cos)
testmn(m,n) = n*nfp
test_fourier = test_fourier .and. &
- check(testmn, result, pol, tor, &
- 'dv_sine_parity_test')
+ & check(testmn, result, pol, tor, &
+ & tor_modes, 'dv_sine_parity_test')
END DO
END DO
@@ -1215,6 +1213,8 @@ FUNCTION test_fourier()
DEALLOCATE(testmn)
DEALLOCATE(result)
+ DEALLOCATE(tor_modes)
+
END FUNCTION
!*******************************************************************************
@@ -1258,27 +1258,32 @@ FUNCTION check_mn(expected, received, m, n, name)
!-------------------------------------------------------------------------------
!> @brief Check all test values.
!>
-!> @param[in] expected Expected value for the test.
-!> @param[in] received Recieved value for the test.
-!> @param[in] name Name of the test.
+!> @param[in] expected Expected value for the test.
+!> @param[in] received Recieved value for the test.
+!> @param[in] mpol Number of poloidal modes.
+!> @param[in] ntor Number of toroidal modes.
+!> @param[in] tor_modes Toroidal modes.
+!> @param[in] name Name of the test.
!-------------------------------------------------------------------------------
- FUNCTION check_all(expected, received, mpol, ntor, name)
+ FUNCTION check_all(expected, received, mpol, ntor, tor_modes, name)
IMPLICIT NONE
! Declare Arguments
- LOGICAL :: check_all
- REAL (rprec), DIMENSION(:,:), INTENT(in) :: expected
- REAL (rprec), DIMENSION(:,:), INTENT(in) :: received
- INTEGER, INTENT(in) :: mpol
- INTEGER, INTENT(in) :: ntor
- CHARACTER (len=*), INTENT(in) :: name
+ LOGICAL :: check_all
+ REAL (rprec), DIMENSION(:,:), INTENT(in) :: expected
+ REAL (rprec), DIMENSION(:,:), INTENT(in) :: received
+ INTEGER, INTENT(in) :: mpol
+ INTEGER, INTENT(in) :: ntor
+ INTEGER, DIMENSION(-ntor:ntor), INTENT(in) :: tor_modes
+ CHARACTER (len=*), INTENT(in) :: name
! Local Variables
- INTEGER :: m
- INTEGER :: n
- INTEGER :: moff
- INTEGER :: noff
+ INTEGER :: m
+ INTEGER :: n
+ INTEGER :: n_mode
+ INTEGER :: moff
+ INTEGER :: noff
! Start of executable code.
moff = LBOUND(expected, 1)
@@ -1286,11 +1291,12 @@ FUNCTION check_all(expected, received, mpol, ntor, name)
check_all = .true.
- DO m = 0, mpol
- DO n = -ntor, ntor
+
+ DO n = -ntor, ntor
+ DO m = 0, mpol
check_all = check_all .and. check(expected(m + moff, n + noff), &
received(m + moff, n + noff), &
- m, n, name)
+ m, tor_modes(n), name)
END DO
END DO
diff --git a/Sources/metrics.f90 b/Sources/metrics.f90
index fc9e4dc..37b1e03 100644
--- a/Sources/metrics.f90
+++ b/Sources/metrics.f90
@@ -122,9 +122,10 @@ END SUBROUTINE init_metric_elements
!>
!> @param[in] mpol_in Number of SIESTA poloidal modes.
!> @param[in] ntor_in Number of SIESTA toroidal modes.
+!> @param[in] nfp_in Number of field periods.
!> @param[in] tor_modes Toroidal mode numbers.
!-------------------------------------------------------------------------------
- SUBROUTINE set_grid_sizes(mpol_in, ntor_in, tor_modes)
+ SUBROUTINE set_grid_sizes(mpol_in, ntor_in, nfp_in, tor_modes)
USE island_params
USE shared_data, ONLY: mblk_size, ndims, lasym
@@ -133,6 +134,7 @@ SUBROUTINE set_grid_sizes(mpol_in, ntor_in, tor_modes)
! Declare Arguments
INTEGER, INTENT(IN) :: mpol_in
INTEGER, INTENT(IN) :: ntor_in
+ INTEGER, INTENT(IN) :: nfp_in
INTEGER, DIMENSION(-ntor_in:ntor_in), INTENT(in) :: tor_modes
! Start of executable code
diff --git a/Sources/perturbation.f90 b/Sources/perturbation.f90
index 4fdfc91..e4470dc 100644
--- a/Sources/perturbation.f90
+++ b/Sources/perturbation.f90
@@ -184,6 +184,7 @@ FUNCTION getwmhd(p)
INTEGER :: iprint
INTEGER :: mres0
INTEGER :: nres0
+ INTEGER :: n
INTEGER :: isign1
INTEGER :: icount
INTEGER :: irscale
@@ -283,13 +284,15 @@ FUNCTION getwmhd(p)
END IF
! Scan in radius to determine primary resonance.
- DO nres0 = -ntor, ntor
+ DO n = -ntor, ntor
! Avoid 0/0 chip/phip at origin
jstart = 3
+ nres0 = fourier_context%tor_modes(n)
+
! Look for multiple resonances
DO WHILE (jstart .lt. ns - 1)
-
+
IF (.not.FindResonance(mres0, nres0, rad, jstart)) THEN
EXIT
END IF
@@ -305,8 +308,9 @@ FUNCTION getwmhd(p)
! vanish at endpoints s=0,1
DO irscale = 1, irad_scale
- CALL GetResPert(irscale, isign1, mres0, nres0, rad, &
- HelPert0, HelPert0A, chip0, phip0, p_width)
+ CALL GetResPert(irscale, isign1, mres0, n, rad, &
+ HelPert0, HelPert0A, chip0, phip0, &
+ p_width)
xc = 0
! FIXME: Why call this again with the same inputs?
wmhd = getwmhd(xc)
@@ -343,7 +347,7 @@ FUNCTION getwmhd(p)
CALL FLUSH(unit_out)
END IF
- CALL GetResPert(irscale, isign1, mres0, nres0, rad, &
+ CALL GetResPert(irscale, isign1, mres0, n, rad, &
HelPert0, HelPert0A, chip0, phip0, p_width)
wmhd = getwmhd(xc)
@@ -503,6 +507,7 @@ SUBROUTINE GetResPert(irscale, isign1, mres0, nres0, rad, &
! Local Variables
INTEGER :: js
INTEGER :: istat
+ INTEGER :: n
REAL (dp) :: rho
REAL (dp) :: rhores
REAL (dp), DIMENSION(:), ALLOCATABLE :: pert_prof
diff --git a/Sources/quantities.f90 b/Sources/quantities.f90
index 2356620..00b48d5 100644
--- a/Sources/quantities.f90
+++ b/Sources/quantities.f90
@@ -485,6 +485,7 @@ SUBROUTINE Init_Bfield(jbsupsmnh, jbsupumnh, jbsupvmnh, &
! local variables
INTEGER :: m
INTEGER :: n
+ INTEGER :: n_mode
INTEGER :: i
INTEGER :: sparity
INTEGER :: mp
@@ -501,7 +502,8 @@ SUBROUTINE Init_Bfield(jbsupsmnh, jbsupumnh, jbsupvmnh, &
DO i = nsmin, nsmax
DO n = -ntor, ntor
- np = n*nfp*sparity
+ n_mode = fourier_context%tor_modes(n)
+ np = n_mode*nfp*sparity
DO m = 0, mpol
mp = m*sparity
jbsupumnh(m,n,i) = phiph(i)*(-np*lmn(m,n,i))
@@ -510,9 +512,9 @@ SUBROUTINE Init_Bfield(jbsupsmnh, jbsupumnh, jbsupvmnh, &
END DO
IF (parity .eq. f_sin) THEN
- jbsupumnh(m0,n0,i) = chiph(i) &
+ jbsupumnh(m0,n0,i) = chiph(i) &
/ fourier_context%orthonorm(m0,n0)
- jbsupvmnh(m0,n0,i) = phiph(i) &
+ jbsupvmnh(m0,n0,i) = phiph(i) &
/ fourier_context%orthonorm(m0,n0)
END IF
END DO
@@ -603,8 +605,10 @@ SUBROUTINE ReCompute_Lambda(lmns, lmnc, jacobh, orthonorm, &
INTEGER :: js
INTEGER :: m
INTEGER :: n
+ INTEGER :: n_mode
INTEGER :: mp
INTEGER :: np
+ INTEGER :: np_mode
INTEGER :: mn
INTEGER :: mnp
INTEGER :: lk
@@ -698,12 +702,15 @@ SUBROUTINE ReCompute_Lambda(lmns, lmnc, jacobh, orthonorm, &
mn = 0
DO n = -ntor, ntor
+ n_mode = fourier_context%tor_modes(n)
DO m = 0, mpol
! Right hand side terms sinmi and -cosmi parts.
mn = mn + 1
- btemp = (-m* (chiph(js)*guv(:,js) + phiph(js)*gvv(:,js)) + &
- n*nfp*(chiph(js)*guu(:,js) + phiph(js)*guv(:,js))) &
+ btemp = (-m* (chiph(js)*guv(:,js) + &
+ & phiph(js)*gvv(:,js)) + &
+ & n_mode*nfp*(chiph(js)*guu(:,js) + &
+ & phiph(js)*guv(:,js))) &
/ jacobh(:,js)
brhs(mn,1) = SUM(cosmni(:,mn)*btemp)
@@ -715,11 +722,14 @@ SUBROUTINE ReCompute_Lambda(lmns, lmnc, jacobh, orthonorm, &
mnp = 0
DO np = -ntor, ntor
+ np_mode = fourier_context%tor_modes(np)
DO mp = 0, mpol
mnp = mnp + 1
! Node the phip term is folded into the solve for lambda.
- atemp = (m*mp*gvv(:,js) + n*nfp*np*nfp*guu(:,js) - &
- nfp*(m*np + mp*n)*guv(:,js))/jacobh(:,js)
+ atemp = (m*mp*gvv(:,js) + &
+ & n_mode*nfp*np_mode*nfp*guu(:,js) - &
+ & nfp*(m*np_mode + mp*n_mode)*guv(:,js)) &
+ & / jacobh(:,js)
amat(mn,mnp) = SUM(cosmni(:,mn)*cosmnp(:,mnp)*atemp)
IF (mnp .eq. mn .and. amat(mn,mnp) .eq. zero) THEN
diff --git a/Sources/siesta_force.f90 b/Sources/siesta_force.f90
index 8a069de..6227400 100644
--- a/Sources/siesta_force.f90
+++ b/Sources/siesta_force.f90
@@ -449,6 +449,7 @@ SUBROUTINE GetMHDForce(fsubsmnf, fsubumnf, fsubvmnf, pijh, &
INTEGER :: sparity
INTEGER :: m
INTEGER :: n
+ INTEGER :: n_mode
INTEGER :: moff
INTEGER :: noff
REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: pmnh
@@ -510,9 +511,11 @@ SUBROUTINE GetMHDForce(fsubsmnf, fsubumnf, fsubvmnf, pijh, &
- m*sparity*pmnf(m,:,nsmin:nsmax)
END DO
DO n = -ntor, ntor
+ n_mode = fourier_context%tor_modes(n)
noff = n + ntor + LBOUND(fsubvmnf, 2)
fsubvmnf(:,noff,nsmin:nsmax) = fsubvmnf(:,noff,nsmin:nsmax) &
- - n*nfp*sparity*pmnf(:,n,nsmin:nsmax)
+ & - (n_mode*nfp*sparity * &
+ & pmnf(:,n,nsmin:nsmax))
END DO
DEALLOCATE(pmnf, pmnf_ds, pmnh)
diff --git a/Sources/siesta_init.f90 b/Sources/siesta_init.f90
index c0f44a4..4167d39 100644
--- a/Sources/siesta_init.f90
+++ b/Sources/siesta_init.f90
@@ -182,11 +182,6 @@ SUBROUTINE init_state(lcurrent_only, lpar_in)
END IF
DEALLOCATE(pfilter)
- CALL CheckFFT(jbsupsmnsh, jbsupumnch, jbsupvmnch, jpmnch, f_sin)
- IF (lasym) THEN
- CALL CheckFFT(jbsupsmnch, jbsupumnsh, jbsupvmnsh, jpmnsh, f_cos)
- END IF
-
IF (.not.lpar) THEN
startglobrow = startr
endglobrow = endr
@@ -498,100 +493,6 @@ SUBROUTINE FilterPressure(parity)
END SUBROUTINE
-!-------------------------------------------------------------------------------
-!> @brief
-!>
-!> @param[in] jbsupsmnh Contravariant magnetic field in the s direction.
-!> @param[in] jbsupumnh Contravariant magnetic field in the u direction.
-!> @param[in] jbsupvmnh Contravariant magnetic field in the v direction.
-!> @param[in] jpmnh Pressure on the half grid.
-!> @param[in] iparity Fourier parity of the s component.
-!-------------------------------------------------------------------------------
- SUBROUTINE CheckFFT(jbsupsmnh, jbsupumnh, jbsupvmnh, jpmnh, iparity)
-
- IMPLICIT NONE
-
-! Declare Arguments
- REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jbsupsmnh
- REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jbsupumnh
- REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jbsupvmnh
- REAL (dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(in) :: jpmnh
- INTEGER, INTENT(in) :: iparity
-
-! Local Variables
- INTEGER :: fours
- INTEGER :: fouruv
- INTEGER :: fcomb
- INTEGER :: sparity
- INTEGER :: m
- INTEGER :: mp
- INTEGER :: n
- INTEGER :: np
- REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: FFT_TEST
- REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: pmnh
- REAL (dp), ALLOCATABLE, DIMENSION(:,:,:) :: temp
- REAL (dp) :: rtest
-
-! Local Parameters
- REAL (dp), PARAMETER :: rtol = 1.E-14_dp
-
-! Start of executable code.
-#undef _TEST_INIT
-!#define _TEST_INIT
-#if defined(_TEST_INIT)
- IF (iparity .EQ. f_sin) THEN
- fours = f_sin
- fouruv = f_cos
- sparity = 1
- ELSE
- fours = f_cos
- fouruv = f_sin
- sparity = -1
- END IF
-
-!check independent variables
- ALLOCATE (fft_test(0:mpol,-ntor:ntor,nsmin:nsmax), &
- temp(ntheta,nzeta,nsmin:nsmax), &
- pmnh(0:mpol,-ntor:ntor,nsmin:nsmax), stat=m)
- CALL assert(m.EQ.0,' Allocation failed in CheckFFT')
- temp = bsupsijh0*jacobh(:,:,nsmin:nsmax)
- CALL tomnsp(temp, fft_test, fours)
- rTest = MAXVAL(ABS(fft_test - jbsupsmnh(:,:,nsmin:nsmax)))
- CALL assert(rtest.LT.rtol, ' jbsupsmnh FFT error in siesta_init')
- temp = bsupuijh0*jacobh(:,:,nsmin:nsmax)
- CALL tomnsp(temp, fft_test, fouruv)
- rTest = MAXVAL(ABS(fft_test - jbsupumnh(:,:,nsmin:nsmax)))
- CALL assert(rTest.LT.rtol,' jbsupumnh FFT error in siesta_init')
- temp = bsupvijh0*jacobh(:,:,nsmin:nsmax)
- CALL tomnsp(temp, fft_test, fouruv)
- rTest = MAXVAL(ABS(fft_test - jbsupvmnh(:,:,nsmin:nsmax)))
- CALL assert(rTest.LT.rtol,' jbsupvmnh FFT error in siesta_init')
- temp = pijh0*jacobh(:,:,nsmin:nsmax)
- CALL tomnsp(temp, fft_test, fouruv)
- rTest = MAXVAL(ABS(fft_test - jpmnh(:,:,nsmin:nsmax)))
- CALL assert(rTest.LT.rtol,' jpmnh FFT error in siesta_init')
-
-! check angle derivatives of pressure
- CALL tomnsp(pijh0, pmnh, fouruv)
- CALL tomnsp(pijh0_du, fft_test, fouruv)
- DO m = 0, mpol
- mp = -sparity*m
- CALL assert(MAXVAL(ABS(fft_test(m,:,:) - mp*pmnh(m,:,:))) .LE. rtol, &
- ' pijh0_du FFT error in siesta_init')
- END DO
-
- CALL tomnsp(pijh0_dv, fft_test, fouruv)
- DO n = -ntor, ntor
- np = -sparity*n*nfp
- CALL assert(MAXVAL(ABS(fft_test(:,n,:) - np*pmnh(:,n,:))) .LE. rtol, &
- ' pijh0_dv FFT error in siesta_init')
- END DO
-
- DEALLOCATE(fft_test, pmnh, temp)
-#endif
-
- END SUBROUTINE
-
!-------------------------------------------------------------------------------
!> @brief Check that p > 0.
!-------------------------------------------------------------------------------
diff --git a/Sources/siesta_run.f90 b/Sources/siesta_run.f90
index db07b2a..60c75dc 100644
--- a/Sources/siesta_run.f90
+++ b/Sources/siesta_run.f90
@@ -92,7 +92,7 @@ FUNCTION siesta_run_construct(run_comm, verbose, init_mpi, close_wout, &
USE siesta_error
USE diagnostics_mod, ONLY: toroidal_flux0
USE blocktridiagonalsolver_s, ONLY: Initialize, GetRanks
- USE siesta_namelist, ONLY: mpolin, ntorin, ntor_modes
+ USE siesta_namelist, ONLY: mpolin, ntorin, nfpin, ntor_modes
USE metrics, ONLY: set_grid_sizes
IMPLICIT NONE
@@ -185,7 +185,7 @@ FUNCTION siesta_run_construct(run_comm, verbose, init_mpi, close_wout, &
siesta_run_construct%time_on = ton
nprecon = init_data(namelist_file)
- CALL set_grid_sizes(mpolin, ntorin, ntor_modes)
+ CALL set_grid_sizes(mpolin, ntorin, nfpin, ntor_modes(-ntorin:ntorin))
CALL second0(toff)
diff --git a/Sources/utilities.f90 b/Sources/utilities.f90
index e5fc42f..27fecfe 100644
--- a/Sources/utilities.f90
+++ b/Sources/utilities.f90
@@ -10,7 +10,7 @@
MODULE utilities
USE stel_kinds
USE island_params, ONLY: ohs=>ohs_i, ns=>ns_i, mpol=>mpol_i, &
- ntor=>ntor_i, nfp=>nfp_i
+ ntor=>ntor_i, nfp=>nfp_i, fourier_context
USE v3_utilities, ONLY: assert, assert_eq
USE fourier, ONLY: m0, m1, m2, n0, n1, f_sin, f_cos, f_ds, f_none, &
b_ds, b_jac, b_con, f_jac
@@ -408,7 +408,7 @@ SUBROUTINE curl_ftoh(asubsmnf, asubumnf, asubvmnf, &
! (sqrt(g)*B^X)mn
DO s = nmin, nsmax
DO n = -ntor, ntor
- np = sparity*n*nfp
+ np = sparity*fourier_context%tor_modes(n)*nfp
DO m = 0, mpol
mp = sparity*m
jbsupsmnh(m,n,s) = np*asubumnh(m,n,s) - mp*asubvmnh(m,n,s)
@@ -551,12 +551,12 @@ SUBROUTINE curl_htof(bsubsmnh, bsubumnh, bsubvmnh, &
CALL assert(nmax .le. UBOUND(jksupsmnf,3), 'UBOUND wrong in curl_htof')
! (sqrt(g)*J^X)mn
- DO m = 0, mpol
- mp = sparity*m
- mo = m + moff
- DO n = -ntor, ntor
- np = sparity*n*nfp
- no = n + noff
+ DO n = -ntor, ntor
+ np = sparity*fourier_context%tor_modes(n)*nfp
+ no = n + noff
+ DO m = 0, mpol
+ mp = sparity*m
+ mo = m + moff
jksupsmnf(mo,no,nsmin:nmax) = np*bsubumnf(m,n,nsmin:nmax) &
- mp*bsubvmnf(m,n,nsmin:nmax)
jksupumnf(mo,no,nsmin:nmax) = np*bsubsmnf(m,n,nsmin:nmax) &
@@ -901,6 +901,7 @@ SUBROUTINE div_h(jbsupsmnh, jbsupumnh, jbsupvmnh, &
INTEGER :: sparity
INTEGER :: m
INTEGER :: n
+ INTEGER :: n_mode
INTEGER :: s
INTEGER :: fours
INTEGER :: fouruv
@@ -930,10 +931,11 @@ SUBROUTINE div_h(jbsupsmnh, jbsupumnh, jbsupvmnh, &
DO s = nsmin, nsmax
DO n = -ntor, ntor
+ n_mode = fourier_context%tor_modes(n)*nfp
DO m = 0, mpol
divmnf(m,n,s) = divmnf(m,n,s) &
- sparity*m*jbsupumnf(m,n,s) &
- - sparity*n*nfp*jbsupvmnf(m,n,s)
+ - sparity*n_mode*jbsupvmnf(m,n,s)
END DO
END DO
END DO