V3FIT
fouri.f
1  SUBROUTINE fouri (grpmn, gsource, amatrix, amatsq, bvec,
2  & bvecNS, ndim)
3  USE vacmod
4  USE parallel_include_module
5  USE timer_sub
6  IMPLICIT NONE
7 C-----------------------------------------------
8 C D u m m y A r g u m e n t s
9 C-----------------------------------------------
10  INTEGER, INTENT(IN) :: ndim
11  REAL(dp), DIMENSION(mnpd2,nuv3), INTENT(IN) :: grpmn
12  REAL(dp), DIMENSION(nuv), INTENT(IN) :: gsource
13  REAL(dp), DIMENSION(mnpd,mnpd,ndim**2), INTENT(OUT) :: amatrix
14  REAL(dp), DIMENSION(mnpd2,mnpd2), INTENT(OUT) :: amatsq
15  REAL(dp), DIMENSION(mnpd,ndim), INTENT(INOUT) :: bvec, bvecNS
16 C-----------------------------------------------
17 C L o c a l P a r a m e t e r s
18 C-----------------------------------------------
19 C interior (int_ext=-1), exterior (int_ext=+1) neumann problem
20  REAL(dp), PARAMETER :: int_ext = 1
21 C-----------------------------------------------
22 C L o c a l V a r i a b l e s
23 C-----------------------------------------------
24  INTEGER :: i, j, k, m, mn, mn0, n
25  REAL(dp), ALLOCATABLE, DIMENSION(:) :: source
26  REAL(dp) :: ton, toff, tfourion, tfourioff
27 !-----------------------------------------------
28 !
29 ! AMATRIX(,1) = A(Sin)(Sin'); AMATRIX(,2) = A(Sin)(Cos');
30 ! AMATRIX(,3) = A(Cos)(Sin'); AMATRIX(,4) = A(Cos)(Cos')
31 !
32 ! ARG OF TRIG FUNCTIONS: m'u' - n'v' CORRESPONDING TO THE PRIMED MESH IN
33 ! PKM EQ.(2.14), AND mu - nv (UNPRIMED MESH) IN EQ. (2.16)
34 !
35 ! ON ENTRY, GRPMN(MN,...) HAS BEEN FOURIER-TRANSFORMED OVER THE UNPRIMED
36 ! COORDINATES (IN FOURP), SO THE FIRST INDEX OF GRPMN CORRESPONDS TO THE FIRST
37 ! INDEX OF AMATRIX. THUS, THE FOURIER TRANSFORMS OF GRPMN HERE ARE OVER THE PRIMED MESH.
38 !
39 ! IN CONTRAST, THE INTEGRAL OF THE SOURCE TERM OVER THE PRIMED MESH WAS ALREADY
40 ! DONE (IN SCALPOT), SO HERE THE FT ARE OVER THE UNPRIMED MESH FOR THE SOURCE.
41 !
42  CALL second0(tfourion)
43 
44  ALLOCATE (source(nuv3), stat=i)
45  IF (i .NE. 0) stop 'Allocation error in fouri'
46 !
47 ! STELLARATOR-SYMMETRIZE SOURCE TERMS (with respect to u,v and -u,-v)
48 ! INDEX (1) IS ANTI-SYMMETRIC, INDEX (2) IS SYMMETRIC
49 !
50 ! GSOURCE = - (2pi)**2 B * dS (h - hsing) * WINT
51 !
52 ! WINT: needed to normalize integral over angles to unity
53 !
54  IF (lasym) THEN
55  source(nuv3min:nuv3max) = onp*gsource(nuv3min:nuv3max)
56  ELSE
57  k = 0
58  DO i = 1, nu2
59  DO j = 1, nv
60  k = k + 1
61  IF (nuv3min.LE.k .AND. k.LE.nuv3max) THEN
62  source(k) = p5*onp*(gsource(k) - gsource(imirr(k)))
63  END IF
64  END DO
65  END DO
66  END IF
67 
68 !
69 ! INITIALIZE RUNNING-SUM ARRAYS
70 !
71  bvecns = 0
72  amatrix = 0
73 
74 !
75 ! PERFORM M,N TRANSFORMS
76 !
77  mn = 0
78  nloop2: DO n = -nf, nf
79  mloop: DO m = 0, mf
80  mn = mn + 1
81  j = 0
82  IF (m.EQ.0 .AND. n.LT.0) cycle
83  DO i = nuv3min, nuv3max
84  j = j + 1
85  bvecns(mn,1) = bvecns(mn,1) + sinmni(j,mn)*source(i)
86 
87  amatrix(:,mn,1) = amatrix(:,mn,1)
88  & + sinmni(j,mn)*grpmn(:mnpd,i)
89 
90  IF (.NOT.lasym) cycle
91 
92  bvecns(mn,2) = bvecns(mn,2) + cosmni(j,mn)*source(i)
93 
94  amatrix(:,mn,2) = amatrix(:,mn,2)
95  & + cosmni(j,mn)*grpmn(:mnpd,i)
96  amatrix(:,mn,3) = amatrix(:,mn,3)
97  & + sinmni(j,mn)*grpmn(mnpd+1:,i)
98  amatrix(:,mn,4) = amatrix(:,mn,4)
99  & + cosmni(j,mn)*grpmn(mnpd+1:,i)
100  END DO
101  END DO mloop
102  END DO nloop2
103 
104  IF (vlactive) THEN
105  CALL second0(ton)
106  CALL mpi_allreduce(mpi_in_place, amatrix, SIZE(amatrix),
107  & mpi_real8, mpi_sum, vac_comm, mpi_err)
108  CALL second0(toff)
109  allreduce_time = allreduce_time + (toff - ton)
110  END IF
111 !
112 ! ADD (still not reduced) ANALYTIC AND NON-SINGULAR PARTS
113 !
114  bvec = bvec + bvecns
115 
116  DEALLOCATE (source, stat=i)
117 
118  mn0 = 1 + mf1*nf !Index of m,n=(0,0)
119 
120 !SANITY CHECKS
121 ! IF (ANY(bvec(1:mn0-mf1:mf1,:) .NE. 0._dp)) STOP 'BVEC != 0'
122 ! IF (ANY(amatrix(:,1:mn0-mf1:mf1,:) .ne. 0._dp)) STOP 'AMAT1 != 0'
123 ! IF (ANY(amatrix(1:mn0-mf1:mf1,:,:) .ne. 0._dp)) STOP 'AMAT2 != 0'
124 !
125 ! ADD DIAGONAL TERMS TO AMATRIX [THE FIRST TERM IN EQ(3.2) OF PKM]
126 !
127  DO mn = 1, mnpd
128  amatrix(mn,mn,1) = amatrix(mn,mn,1) + pi3*int_ext
129  END DO
130 
131  IF (lasym) THEN
132  DO mn = 1, mnpd
133  amatrix(mn,mn,4) = amatrix(mn,mn,4) + pi3*int_ext
134  END DO
135  amatrix(mn0,mn0,4) = amatrix(mn0,mn0,4) + pi3*int_ext !!COS(0-0) mode *2
136  END IF
137 
138 !
139 ! PUT ELEMENTS INTO SQUARE MATRIX
140 !
141  amatsq(:mnpd,:mnpd) = amatrix(:,:,1) !Sin-Sin'
142 
143  IF (lasym) THEN
144  amatsq(:mnpd,1+mnpd:mnpd2) = amatrix(:,:,2) !Sin-Cos'
145  amatsq(1+mnpd:mnpd2,:mnpd) = amatrix(:,:,3) !Cos-Sin'
146  amatsq(1+mnpd:mnpd2,1+mnpd:mnpd2) = amatrix(:,:,4) !Cos-Cos'
147  END IF
148 
149  CALL second0(tfourioff)
150  timer_vac(tfouri) = timer_vac(tfouri) + (tfourioff-tfourion)
151  fouri_time = timer_vac(tfouri)
152 
153  END SUBROUTINE fouri