V3FIT
analyt.f
1  SUBROUTINE analyt (grpmn, bvec, ivacskip, ndim)
2  USE vacmod
3  USE parallel_include_module
4  USE timer_sub
5  IMPLICIT NONE
6 C-----------------------------------------------
7 C D u m m y A r g u m e n t s
8 C-----------------------------------------------
9  INTEGER, INTENT(IN) :: ivacskip, ndim
10  REAL(dp), INTENT(OUT) :: grpmn(mnpd2,nuv3)
11  REAL(dp), INTENT(OUT) :: bvec(mnpd,ndim)
12 C-----------------------------------------------
13 C L o c a l V a r i a b l e s
14 C-----------------------------------------------
15  INTEGER :: l, n, m, i, q, j, k, ll, blksize, mn
16  REAL(dp), DIMENSION(:), ALLOCATABLE ::
17  1 r0p, r1p, r0m, r1m, sqrtc, sqrta, tlp2, tlp1, tlp, tlm2,
18  2 tlm1, tlm, adp, adm, cma, ra1p, ra1m, slm, slp, tlpm, slpm,
19  3 delt1u, azp1u, azm1u, cma11u, sqad1u, sqad2u
20  REAL(dp) :: fl, fl1, fl2, sign1, tanalon, tanaloff
21 C-----------------------------------------------
22  CALL second0(tanalon)
23 
24  ALLOCATE (r0p(nuv3), r1p(nuv3), r0m(nuv3), r1m(nuv3),
25  1 sqrtc(nuv3), sqrta(nuv3), tlp2(nuv3), tlp1(nuv3),
26  2 tlp(nuv3), tlm2(nuv3), tlm1(nuv3), tlm(nuv3), adp(nuv3),
27  3 adm(nuv3), cma(nuv3), ra1p(nuv3), ra1m(nuv3), slm(nuv3),
28  4 slp(nuv3), tlpm(nuv3), slpm(nuv3), delt1u(nuv3),
29  5 azp1u(nuv3), azm1u(nuv3), cma11u(nuv3), sqad1u(nuv3),
30  6 sqad2u(nuv3), stat = l)
31  IF (l .ne. 0) stop 'Allocation error in SUBROUTINE analyt'
32 
33 !
34 ! ALL EQUATIONS REFER TO THE PAPER BY P. MERKEL (PKM)
35 ! IN J. COMPUT. PHYSICS 66, p83 (1986)
36 !
37 ! IN GOING BETWEEN THE COMPLEX NOTATION OF (PKM) AND OUR REAL FORM,
38 ! NOTE THAT THE INTEGRALS (APPENDIX, PKM) Imn AND Kmn ARE BOTH REAL.
39 ! THUS, THE SIN(mu-nv) INTEGRALS OF THE SINGULAR PIECE (ANALYTIC CONTRIBUTION)
40 ! VANISHES.
41 !
42 ! THE REQUIRED SOURCE-TERM INTEGRALS ARE (Eq.2.16-2.17):
43 !
44 ! BVECS(m,n) = Int< SIN(mu' - nv') han(u',v') >
45 ! BVECC(m,n) = Int< COS(mu' - nv') han(u',v') >
46 !
47 ! Where Int<...> means integration over u (theta) and v (zeta) and
48 ! summation over field periods. These can be written in terms of PKM integrals
49 ! Imn(a,b,c), where a(u,v) = guu (g theta-theta), etc.:
50 !
51 ! BVECS(m,n) = ALP * Int<SIN(mu' - nv') * F * Im,-n(a,b,c)>
52 ! BVECC(m,n) = ALP * Int<COS(mu' - nv') * F * Im,-n(a,b,c)>
53 !
54 ! Here, F = - BNORM(u',v') is defined in Eq.(2.13), and ALP = (2*pi/nfp).
55 !
56 ! Similarly, the analytic part of the matrix A(m,n;m',n') can be written:
57 !
58 ! A(m,n;m',n') = (2*pi/nfp) * Int<SIN(mu' - nv')*SIN(m'u' - n'v')
59 ! [Km,-n](a',b',c';A',B',C')>
60 !
61 ! On EXIT, GRPMN(ip,m,n) = ALP * SIN(ip,m,n) * K[m,-n](ip)
62 !
63 !
64 ! COMPUTE ALL QUANTITIES INDEPENDENT OF THE MODE INDICES L,M,N
65 ! NOTE: 2b = guv_b HAS FACTOR OF 2 BUILT IN (see SUBROUTINE SURFACE)
66 !
67 ! ADP(M): a +(-)2b + c
68 ! CMA: c - a
69 ! DELTA: 4*(ac - b**2)
70 ! AZP(M): A +(-)2*B + C
71 ! CMA1: C - A
72 ! R1P(M): Coefficient of l*Tl+(-) in eq (A17)
73 ! R0P(M): Coefficient of l*T(l-1)+(-) in eq (A17)
74 ! RA1P(M):Coefficient of Tl+(-) in eq (A17)
75 !
76  DO k = nuv3min, nuv3max
77  adp(k) = guu_b(k) + guv_b(k) + gvv_b(k)
78  adm(k) = guu_b(k) - guv_b(k) + gvv_b(k)
79  cma(k) = gvv_b(k) - guu_b(k)
80  sqrtc(k) = two*sqrt(gvv_b(k))
81  sqrta(k) = two*sqrt(guu_b(k))
82  END DO
83 
84  IF (ivacskip .EQ. 0) THEN
85 
86  grpmn(:,nuv3min:nuv3max) = 0
87 
88  DO k = nuv3min, nuv3max
89  delt1u(k) = adp(k)*adm(k) - cma(k)*cma(k)
90  azp1u(k) = auu(k) + auv(k) + avv(k)
91  azm1u(k) = auu(k) - auv(k) + avv(k)
92  cma11u(k)= avv(k) - auu(k)
93  r1p(k) = (azp1u(k)*(delt1u(k) - cma(k)*cma(k))/adp(k)
94  1 - azm1u(k)*adp(k) + two*cma11u(k)*cma(k))/delt1u(k)
95  r1m(k) = (azm1u(k)* (delt1u(k) - cma(k)*cma(k))/adm(k)
96  1 - azp1u(k)*adm(k) + two*cma11u(k)*cma(k))/delt1u(k)
97  r0p(k) = (-azp1u(k)*adm(k)*cma(k)/adp(k) - azm1u(k)*cma(k)
98  1 + two*cma11u(k)*adm(k))/delt1u(k)
99  r0m(k) = (-azm1u(k)*adp(k)*cma(k)/adm(k) - azp1u(k)*cma(k)
100  1 + two*cma11u(k)*adp(k))/delt1u(k)
101  ra1p(k) = azp1u(k)/adp(k)
102  ra1m(k) = azm1u(k)/adm(k)
103  END DO
104  ENDIF
105 
106 !
107 ! INITIALIZE VECTORS
108 !
109 ! bvec = 0
110 !
111 ! INITIALIZE T0+ and T0-
112 !
113 ! TLP(M): TL+(-)
114 ! TLP(M)1:T(L-1)+(-)
115 ! TLP(M)2:T(L-2)+(-)
116 !
117  DO k = nuv3min,nuv3max
118  sqad1u(k) = sqrt(adp(k))
119  sqad2u(k) = sqrt(adm(k))
120  tlp1(k) = 0
121  tlm1(k) = 0
122  tlp(k) = one/sqad1u(k)*log((sqad1u(k)*sqrtc(k)
123  1 + adp(k) + cma(k))/(sqad1u(k)*sqrta(k)
124  2 - adp(k) + cma(k)))
125  tlm(k) = one/sqad2u(k)*log((sqad2u(k)*sqrtc(k)
126  1 + adm(k) + cma(k))/(sqad2u(k)*sqrta(k)
127  2 - adm(k) + cma(k)))
128  tlpm(k) = tlp(k) + tlm(k)
129  END DO
130 !
131 ! BEGIN L-SUM IN EQ (A14) TO COMPUTE Imn (and Kmn) INTEGRALS
132 ! NOTE THAT IN THE LOOP OVER L BELOW: L == |m - n| + 2L_A14
133 ! THUS, L BELOW IS THE INDEX OF THE T+- (S+-)
134 !
135  sign1 = 1
136  fl1 = 0
137 
138  lloop: DO l = 0, mf + nf
139  fl = fl1
140 !
141 ! COMPUTE SL+ and SL- , Eq (A17)
142 ! SLP(M): SL+(-)
143 !
144  IF (ivacskip .eq. 0) THEN
145  DO k = nuv3min,nuv3max
146  slp(k) = (r1p(k)*fl + ra1p(k))*tlp(k) + r0p(k)*fl*tlp1(k)
147  1 - (r1p(k) + r0p(k))/sqrtc(k)
148  2 + sign1*(r0p(k) - r1p(k))/sqrta(k)
149  slm(k) = (r1m(k)*fl + ra1m(k))*tlm(k) + r0m(k)*fl*tlm1(k)
150  1 - (r1m(k) + r0m(k))/sqrtc(k)
151  2 + sign1*(r0m(k) - r1m(k))/sqrta(k)
152  slpm(k) = slp(k) + slm(k)
153  END DO
154  ENDIF
155 !
156 ! BEGIN MODE NUMBER (m,n) LOOP
157 !
158  DO n = 0, nf
159  DO m = 0, mf
160 
161  IF (l .EQ. 0) THEN
162  mn = m + mf1*(n+nf) + 1
163  bvec(mn,:) = 0
164  mn = m + mf1*(nf-n) + 1
165  bvec(mn,:) = 0
166  END IF
167 
168  IF (cmns(l,m,n) .eq. zero) cycle
169 
170  IF (n.eq.0 .or. m.eq.0) THEN
171 !
172 ! 1. n = 0 and m >= 0 OR n > 0 and m = 0
173 !
174  CALL analysum (grpmn, bvec, slpm, tlpm, m, n, l,
175  1 ivacskip, ndim)
176 
177  ELSE
178 !
179 ! 2. n>=1 and m>=1
180 !
181  CALL analysum2 (grpmn, bvec, slm, tlm, slp, tlp,
182  1 m, n, l, ivacskip, ndim)
183 
184  ENDIF
185  END DO
186  END DO
187 
188 !
189 ! UPDATE "TL's" (FOR L -> L+1) USING EQ (A15)
190 !
191  fl1 = fl1 + 1 !next l
192  fl2 = 2*fl1-1
193  sign1 = -sign1 !(-1)**l (next l now)
194  DO k = nuv3min, nuv3max
195  tlp2(k) = tlp1(k)
196  tlm2(k) = tlm1(k)
197  tlp1(k) = tlp(k)
198  tlm1(k) = tlm(k)
199  tlp(k) = ((sqrtc(k) + sign1*sqrta(k)) - fl2*
200  1 cma(k)*tlp1(k) - fl*adm(k)*tlp2(k))/(adp(k)*fl1)
201  tlm(k) = ((sqrtc(k) + sign1*sqrta(k)) - fl2*
202  1 cma(k)*tlm1(k) - fl*adp(k)*tlm2(k))/(adm(k)*fl1)
203  tlpm(k) = tlp(k) + tlm(k)
204  END DO
205 
206  END DO lloop
207 
208  DEALLOCATE (r0p, r1p, r0m, r1m, sqrtc, sqrta, tlp2, tlp1,
209  1 tlp, tlm2, tlm1, tlm, adp, adm, cma, ra1p, ra1m, slm,
210  2 slp, tlpm, slpm, delt1u, azp1u, azm1u, cma11u, sqad1u,
211  3 sqad2u, stat = l)
212 
213  CALL second0(tanaloff)
214  timer_vac(tanal) = timer_vac(tanal) + (tanaloff-tanalon)
215  analyt_time = timer_vac(tanal)
216 
217 
218  END SUBROUTINE analyt
219