V3FIT
getbrho.f
1  SUBROUTINE getbsubs (bsubsmn, frho, bsupu, bsupv, mmax,
2  1 nmax, info)
3  USE stel_kinds
4  USE vmec_input, ONLY: nfp, nzeta, lasym
5  USE vmec_dim, ONLY: ntheta1, ntheta2, ntheta3
6  USE vmec_persistent, ONLY: cosmu, sinmu, cosnv, sinnv
7  USE vmec_main, ONLY: r0scale
8  IMPLICIT NONE
9 C-----------------------------------------------
10 C D u m m y A r g u m e n t s
11 C-----------------------------------------------
12  INTEGER, INTENT(in) :: mmax, nmax
13  INTEGER, INTENT(out) :: info
14  REAL(rprec), INTENT(out) :: bsubsmn(0:mmax, -nmax:nmax, 0:1)
15  REAL(rprec), DIMENSION(nzeta, ntheta3), INTENT(in) ::
16  1 bsupu, bsupv, frho
17 C-----------------------------------------------
18 C L o c a l P a r a m e t e r s
19 C-----------------------------------------------
20  REAL(rprec), PARAMETER :: p5 = 0.5_dp, one = 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, m, n, nmax1, itotal, ijtot, mntot
25  REAL(rprec) :: ccmn, ssmn, csmn, scmn, dm, dn, termsc, termcs,
26  ` termcc, termss, amn
27  REAL(rprec), ALLOCATABLE :: amatrix(:,:), save_matrix(:,:),
28  1 brhs(:)
29  LOGICAL :: lpior0
30  EXTERNAL solver
31 C-----------------------------------------------
32 !
33 ! Solves the radial force balance B dot bsubs = Fs for bsubs in real space using collocation
34 ! Here, Fs = frho(mn) is the Fourier transform SQRT(g)*F (part of radial force
35 ! balance sans the B dot bsubs term)
36 !
37 ! Storage layout: bsubsmn(0:mmax, 0:nmax,0) :: coefficient of sin(mu)cos(nv)
38 ! bsubsmn(0:mmax,-1:nmax,0) :: coefficient of cos(mu)sin(nv)
39 ! for lasym = T
40 ! bsubsmn(0:mmax, 0:nmax,1) :: coefficient of cos(mu)cos(nv)
41 ! bsubsmn(0:mmax,-1:nmax,1) :: coefficient of sin(mu)sin(nv)
42 !
43 ! where 0<=m<=mmax and 0<=n<=nmax
44 !
45  info = -3
46  IF ((mmax+1 .ne. ntheta2) .or. (nmax .ne. nzeta/2)) RETURN
47 
48  nmax1 = max(0,nmax-1)
49  itotal = ntheta3*nzeta
50  IF (.not.lasym) itotal = itotal - 2*nmax1
51  ALLOCATE (amatrix(itotal, itotal),
52  1 brhs(itotal), save_matrix(itotal, itotal), stat=m)
53  IF (m .ne. 0) stop 'Allocation error in getbsubs'
54 
55  amatrix = 0
56 
57 !
58 ! bsubs = BSC(M,N)*SIN(MU)COS(NV) + BCS(M,N)*COS(MU)SIN(NV)
59 ! + BCC(M,N)*COS(MU)COS(NV) + BSS(M,N)*SIN(MU)SIN(NV) (LASYM=T ONLY)
60 !
61 
62  ijtot = 0
63  brhs = 0
64 
65  DO i = 1, ntheta3
66  DO j = 1, nzeta
67 ! IGNORE u=0,pi POINTS FOR v > pi: REFLECTIONAL SYMMETRY
68  lpior0 = ((i.eq.1 .or. i.eq.ntheta2) .and. (j.gt.nzeta/2+1))
69  IF (lpior0 .and. .not. lasym) cycle
70  ijtot = ijtot + 1
71  brhs(ijtot) = frho(j,i)
72  mntot = 0
73  DO m = 0, mmax
74  DO n = 0, nmax
75  IF (mntot .ge. itotal) EXIT
76  IF (m.eq.0 .and. n.eq.0 .and. lasym) cycle
77  mntot = mntot+1
78  ccmn = cosmu(i,m)*cosnv(j,n)
79  ssmn = sinmu(i,m)*sinnv(j,n)
80  dm = m * bsupu(j,i)
81  dn = n * bsupv(j,i) * nfp
82  termsc = dm*ccmn - dn*ssmn
83  termcs =-dm*ssmn + dn*ccmn
84  IF (n.eq.0 .or. n.eq.nmax) THEN
85  IF (m .gt. 0) THEN
86  amatrix(ijtot,mntot) = termsc !!ONLY bsc != 0 for n=0, nmax1
87  ELSE IF (n .eq. 0) THEN
88  amatrix(ijtot,mntot) = bsupv(j,i) !!pedestal for m=0,n=0 mode, which should = 0
89  ELSE
90  amatrix(ijtot,mntot) = termcs !!bcs(m=0,n=nmax)
91  END IF
92  ELSE IF (m.eq.0 .or. m.eq.mmax) THEN
93  amatrix(ijtot,mntot) = termcs !!ONLY bcs != 0 for m=0,mmax
94  ELSE
95  amatrix(ijtot,mntot) = termsc
96  mntot = mntot+1
97  amatrix(ijtot,mntot) = termcs
98  END IF
99 
100  IF (.not.lasym) cycle
101  IF (m.eq.0 .and. (n.eq.0 .or. n.eq.nmax)) cycle
102 
103  IF (mntot .ge. itotal) EXIT
104  mntot = mntot+1
105  csmn = cosmu(i,m)*sinnv(j,n)
106  scmn = sinmu(i,m)*cosnv(j,n)
107  termcc =-dm*scmn - dn*csmn
108  termss = dm*csmn + dn*scmn
109 
110  IF ((n.eq.0 .or. n.eq.nmax) .or.
111  1 (m.eq.0 .or. m.eq.mmax)) THEN
112  amatrix(ijtot,mntot) = termcc !!ONLY bcc != 0 for m=0 or mmax
113  ELSE
114  amatrix(ijtot,mntot) = termcc
115  mntot = mntot+1
116  amatrix(ijtot,mntot) = termss
117  END IF
118 
119  END DO
120  END DO
121  END DO
122  END DO
123 
124  save_matrix = amatrix
125 
126  info = -1
127  IF (ijtot .ne. itotal .or. mntot .ne. itotal) THEN
128  print *,' itotal = ', itotal,' ijtot = ', ijtot,
129  1 ' mntot = ', mntot
130  print *,' ntheta3: ',ntheta3,' nzeta: ', nzeta,
131  1 ' mnyq: ', mmax,' nnyq: ', nmax
132  GOTO 200
133  ELSE
134  CALL solver (amatrix, brhs, itotal, 1, info)
135  IF (info .ne. 0) GOTO 200
136  END IF
137 
138 !
139 ! CHECK SOLUTION FROM SOLVER
140 !
141 ! GOTO 100
142  ijtot = 0
143  DO i = 1, ntheta3
144  DO j = 1, nzeta
145  lpior0 = ((i.eq.1 .or. i.eq.ntheta2) .and. (j.gt.nzeta/2+1))
146  IF (lpior0 .and. .not.lasym) cycle
147  ijtot = ijtot + 1
148  amn = sum(save_matrix(ijtot,:)*brhs(:))
149  IF (abs(amn) .lt. 1.e-12_dp) cycle
150  IF (abs(frho(j,i) - amn) .gt. 1.e-8_dp*abs(amn)) THEN
151  print 50,'In GETbsubs, i = ',i,' j = ',j,
152  1 ' Original force = ', frho(j,i),' Final force = ', amn
153  END IF
154  END DO
155  END DO
156  50 FORMAT(a,i5,a,i5,a,1p,e10.3,a,1p,e10.3)
157  100 CONTINUE
158 !
159 ! CONVERT BACK TO BS*SIN(MU - NV) REPRESENTATION
160 ! AND (FOR lasym) BC*COS(MU - NV)
161 !
162  mntot = 0
163  bsubsmn = 0
164  DO m = 0, mmax
165  DO n = 0, nmax
166  IF (mntot .ge. itotal) EXIT
167  IF (m.eq.0 .and. n.eq.0 .and. lasym) cycle
168  mntot = mntot+1
169  IF (n.eq.0 .or. n.eq.nmax) THEN
170  IF (m .gt. 0) THEN
171  bsubsmn(m,n,0) = brhs(mntot)
172  ELSE IF (n .eq. 0) THEN
173  bsubsmn(m,n,0) = brhs(mntot)
174  ELSE
175  bsubsmn(m,-n,0) = brhs(mntot)
176  END IF
177  ELSE IF (m.eq.0 .or. m.eq.mmax) THEN
178  bsubsmn(m,-n,0) = brhs(mntot)
179  ELSE
180  bsubsmn(m,n,0) = brhs(mntot)
181  mntot = mntot+1
182  bsubsmn(m,-n,0) = brhs(mntot)
183  END IF
184 
185  IF (.not.lasym) cycle
186  IF (m.eq.0 .and. (n.eq.0 .or. n.eq.nmax)) cycle
187  IF (mntot .ge. itotal) EXIT
188  mntot = mntot+1
189 
190  IF ((n.eq.0 .or. n.eq.nmax) .or.
191  1 (m.eq.0 .or. m.eq.mmax)) THEN
192  bsubsmn(m,n,1) = brhs(mntot)
193  ELSE
194  bsubsmn(m,n,1) = brhs(mntot)
195  mntot = mntot+1
196  bsubsmn(m,-n,1)= brhs(mntot)
197  END IF
198 
199  END DO
200  END DO
201 
202  IF (mntot .ne. ijtot) info = -2
203 
204  200 CONTINUE
205 
206  DEALLOCATE (amatrix, save_matrix, brhs)
207 
208  END SUBROUTINE getbsubs
Generated on Thu Mar 5 2020 15:49:27 for V3FIT by   doxygen 1.8.17