V3FIT
descur_sub.f
1  SUBROUTINE descur_sub(g1,g2)
2 c THIS IS PROGRAM DESCUR - WHICH USES A STEEPEST DESCENT
3 c ALGORITHM TO FIND A LEAST SQUARES APPROXIMATION TO AN
4 c ARBITRARY 3-D SPACE CURVE. ANGLE CONSTRAINTS, BASED ON A
5 c MINIMUM SPECTRAL WIDTH CRITERION, ARE APPLIED.
6 c THE CONSTRAINT IS SATISFIED BY TANGENTIAL VARIATIONS
7 c ALONG THE CURVE.
8 c
9 c THE MAIN PROGRAM SETS UP THE INITIAL POINTS TO BE FITS AND
10 c THEN CALLS THE SUBROUTINE SCRUNCH, WHICH DOES THE ACTUAL
11 c CURVE-FITTING.
12 c
13 c***********************************************************************
14 c REVISION 1: January 26, 1989
15 c Compute Fit to individual toroidal planes separately
16 c and THEN Fourier decompose in phi, using optimized
17 c angle representation
18 c
19 c REVISION 2: July, 1995
20 c Generalized for up-down asymmetric CASE
21 c
22 c REVISION 3: July, 1997
23 c Converted to Fortran-90
24 c Eliminated Lagrange multiplier constraint
25 c by using an EXPlicit representation
26 c
27 c
28 c PARAMETER VALUE INTERPRETATIONS:
29 c
30 c MPOL: ACTUAL NO. OF POLOIDAL MODES USED IN FIT FOR R and Z
31 c MRHO: ACTUAL NO. OF POLOIDAL MODES IN FIT FOR QUASI-POLAR RHO
32 c NTHETA: ACTUAL NO. OF THETA (POLOIDAL) POINTS AUSED IN FIT
33 c NPHI: ACTUAL NO. OF TOROIDAL PLANES IN FIT
34 c
35 c***********************************************************************
36 C-----------------------------------------------
37 C M o d u l e s
38 C-----------------------------------------------
39  USE vname0
40  USE vname1
41  USE gfile
42  USE mapout
43  IMPLICIT NONE
44 C-----------------------------------------------
45 C L o c a l P a r a m e t e r s
46 C-----------------------------------------------
47  INTEGER, PARAMETER :: nphi20 = 1 + nv/2
48  INTEGER, PARAMETER :: mntot = mu*(2*nphi20 + 1)
49 C-----------------------------------------------
50 C L o c a l V a r i a b l e s
51 C-----------------------------------------------
52  INTEGER :: niter, nstep, itype,
53  1 i, k, j, m, n, mexp
54  REAL(rprec), DIMENSION(nu*nv) :: rin, zin, LENgth
55  REAL(rprec), DIMENSION(mntot) :: rbc, rbs, zbc, zbs
56  REAL(rprec), DIMENSION(0:nv) :: rmnaxis, zmnaxis
57  REAL(rprec), DIMENSION(0:mu-1,-nphi20:nphi20) ::
58  1 rbdy_3d, zbdy_3d
59  REAL(rprec), DIMENSION(4) :: rbdy, zbdy, rdshape, zdshape,
60  1 rbean, zbean, rbelt, zbelt, rellip, zellip, rsqr, zsqr
61  REAL(rprec), DIMENSION(0:3,0:2) :: rheliac, zheliac
62  REAL(rprec) :: ftol, pexp, qexp, zeta, texp,
63  1 dumr, dumphi, dumz, arg, PSI, XI, DELTA_R,
64  2 rmin, rmax, denom
65  CHARACTER :: datafile*80
66  CHARACTER*1 :: ch_yn
67 C-----------------------------------------------
68 
69  TYPE(GEQDSK) :: g1
70  TYPE(RSZS) :: g2
71  data rdshape/3.0, 0.991, 0.136, 0./
72  data zdshape/0.0, 1.409, -0.118, 0./
73  data rbean/3.0, 1.042, 0.502, - 0.0389/
74  data zbean/0.0, 1.339, 0.296, - 0.0176/
75  data rbelt/3.0, 0.453, 0., 0./
76  data zbelt/0.0, 0.60, 0., 0.196/
77  data rellip/3.0e0, 1.000e0, 0., 0./
78  data zellip/0.0, 6.00e0, 0., 0./
79  data rheliac/4.115, 0.4753, 0., 0., 0.3225, -0.06208, -0.1358,
80  1 -0.04146, 0.0136, 0.0806, -0.0205, 0.0445/
81  data zheliac/0., -0.5045, 0., 0., 0.337, 0.0637, 0.1325, -0.04094
82  1 , 0.01152, 0.06186, -0.03819, - 0.02366/
83  data rsqr/ 3.0, 0.4268, 0., 0.07322/,
84  1 zsqr/ 0.0, 0.4268, 0., -0.07322/
85 
86 c**********************************************************************
87 c CONTROL DATA THAT CAN BE SPECIFIED BY USER
88 c**********************************************************************
89  data niter/1500/
90  data nstep/100/
91 
92  data ftol/1.e-5_dbl/
93  data pexp/4.0_dbl/
94  data qexp/1.0_dbl/
95  data mexp/4/
96 
97  data datafile/'none'/
98 
99  twopi = 8*atan(one)
100 
101 c**********************************************************************
102 c CREATE ARRAYS OF CURVE POINTS, RIN(THETA,PHI) , ZIN(THETA,PHI),
103 c TO BE FIT TO AN OPTIMIZED FOURIER SERIES BY CALLING "SCRUNCH".
104 c NOTE: N IS IN UNITS NORMED TO NFP (NUMBER OF FIELD PERIODS) IN 20 LOOP
105 c**********************************************************************
106 c
107 c Read in data
108 c
109 
110  hb_parameter = 1
111 
112  nfp = 1
113 
114  itype = 0
115 
116  1002 CONTINUE
117  ntheta = g2%nthet ; nphi = 1 ; nfp = 1
118  IF (ntheta .gt. nu) stop 'ntheta > nu'
119  IF (nphi .gt. nv) stop 'nphi > nv'
120  mpol = mu
121  mpol1 = mpol - 1
122  nphi2 = 1 + nphi/2
123  i = 0 ; rin = 0 ; zin = 0
124 ! rin(1:ntheta)=g2%rs(1:ntheta,1)
125 ! zin(1:ntheta)=g2%zs(1:ntheta,1)
126  rin(1:ntheta)=g2%rs(g2%npsi,1:ntheta)
127  zin(1:ntheta)=g2%zs(g2%npsi,1:ntheta)
128  4000 CONTINUE
129 c**********************************************************************
130 c PERFORM OPTIMIZING CURVE-FIT
131 c**********************************************************************
132  CALL scrunch (rin, zin, pexp, qexp, rbc, zbs, rbs, zbc, rmnaxis,
133  1 zmnaxis, ftol, niter, nstep, mexp)
134 
135 c**********************************************************************
136 c PRINT OPTIMIZED TRANSFORM COEFFICIENTS TO "OUTCURVE" FILE
137 c AND STORE RIN, ZIN, RBC, ZBS DATA IN PLOTOUT FILE
138 c**********************************************************************
139  CALL printit(rin,zin,rbc,zbs,rbs,zbc,rmnaxis,zmnaxis,g2)
140  CLOSE(10)
141 
142  END SUBROUTINE descur_sub
143  SUBROUTINE scrunch(rin, zin, pexp, qexp, rbc, zbs, rbs, zbc,
144  1 rmnaxis, zmnaxis, ftol, niter, nstep, mexp)
145 C-----------------------------------------------
146 C M o d u l e s
147 C-----------------------------------------------
148  USE vname0
149  USE vname1
150  IMPLICIT NONE
151 C-----------------------------------------------
152 C D u m m y A r g u m e n t s
153 C-----------------------------------------------
154  INTEGER :: niter, nstep, mexp
155  REAL(rprec) pexp, qexp, ftol
156  REAL(rprec), DIMENSION(ntheta,nphi) :: rin, zin
157  REAL(rprec), DIMENSION(*) :: rbc, zbs, rbs, zbc,
158  1 rmnaxis, zmnaxis
159 C-----------------------------------------------
160 C L o c a l V a r i a b l e s
161 C-----------------------------------------------
162  INTEGER :: irst=1, nplane, imodes, iter, modeno, mrho1
163  REAL(rprec), ALLOCATABLE, DIMENSION(:,:) :: angle
164  REAL(rprec), DIMENSION(2*mpol,nphi) :: result
165  REAL(rprec) :: too_large,
166  1 gtrig, fsq, gmin, g11, gout
167 C-----------------------------------------------
168 c
169 c PARAMETER VALUES TO BE SET IN NAME1.INC:
170 c
171 c MU: MAXIMUM NO. OF POLOIDAL MODES IN FIT
172 c NU: MAXIMUM NO. OF THETA (POLOIDAL) POINTS ALLOWED IN FIT
173 c NV: MAXIMUM NO. OF TOROIDAL PLANES IN FIT
174 c
175 !
176 
177 !
178 ! CHECK THAT INPUT DIMENSIONS ARE CONSISTENT WITH NAME1 PARAMETERS
179 !
180  deltf = 0.97_dbl
181  IF (mod(nphi,2).ne.0 .and. nphi.ne.1 .or. nphi.eq.0) THEN
182  print *, ' NPHI > 0 must be EVEN for non-symmetric systems'
183  stop
184  ENDIF
185  IF (ntheta > nu) THEN
186  print *, ' NTHETA input to SCRUNCH is too large'
187  stop
188  ENDIF
189  IF (nphi > nv) THEN
190  print *, ' NPHI input to SCRUNCH is too large'
191  stop
192  ENDIF
193  IF (mpol > mu) THEN
194  print *, ' MPOL input to SCRUNCH is too large'
195  stop
196  ENDIF
197  IF (nfp .le. 0) THEN
198  print *, ' NFP > 0 must be positive and exceed zero'
199  stop
200  ENDIF
201 
202 
203 !
204 ! INITIALIZE FIXED m,n ARRAYS
205 !
206 CEAL OPEN(unit=3, file='outcurve', status='unknown')
207  twopi = 8*atan(one)
208  CALL fixaray (pexp, qexp, mexp)
209 
210 !
211 ! COMPUTE INITIAL GUESSES (MUST ORDER ANGLES MONOTONICALLY FOR
212 ! NON-STARLIKE DOMAINS)
213 !
214  ALLOCATE (angle(ntheta, nphi))
215  CALL inguess (rin, zin, angle)
216 
217  too_large = 1.e6_dbl
218  gtrig = 1.e-4_dbl
219  mrho1 = mrho+1
220 
221  ALLOCATE (xvec(n2), gvec(n2), xdot(n2), xstore(n2))
222 !
223 ! BEGIN MAIN INTERATION LOOP
224 !
225  DO nplane = 1, nphi
226 !
227 ! INITIALIZE M = 0 and M = 1 MODE AMPLITUDES
228 !
229 ! STACKING OF XVEC ARRAY (FOR FIXED TOROIDAL PLANE)
230 ! XVEC(1) : R0c (Symmetric components)
231 ! XVEC(2,mrho1) : Rhoc (Symmetric components)
232 ! XVEC(1+mrho1) : Z0c (Asymmetric components)
233 ! XVEC(2+mrho1,2*mrho1): Rhos (Asymmetric components)
234 ! XVEC(2*mrho1+1,n2) : Theta angle
235 !
236  xstore(:n2) = zero
237  xdot(:n2) = zero
238  xvec(:n2) = zero
239  r0n = zero ; z0n = zero ; raxis = zero ; zaxis = zero
240 CEAL WRITE (3, 10) nplane
241 CEAL WRITE (*, 10) nplane
242  10 FORMAT(/' Fitting toroidal plane # ',i3)
243 
244  CALL amplitud (r0n(nplane), z0n(nplane), angle(1,nplane),
245  1 xvec(1), xvec(1+mrho1), xvec(2:mrho1), xvec(2+mrho1:2*mrho1),
246  2 xvec(2*mrho1+1:), rin(1,nplane), zin(1,nplane))
247 
248 CEAL WRITE (3, 20)
249 CEAL WRITE (*, 20)
250  20 FORMAT(/' ITERATIONS RMS ERROR FORCE GRADIENT <M>',
251  1 ' MAX m DELT')
252 
253  imodes = mrho
254  delt = deltf
255 
256  DO iter = 1, niter
257  gvec(:n2) = zero
258  CALL funct (xvec(1), xvec(1+mrho1), xvec(2:mrho1),
259  1 xvec(2+mrho1:2*mrho1),xvec(2*mrho1+1:), gvec(1),
260  2 gvec(1+mrho1), gvec(2:mrho1), gvec(2+mrho1:2*mrho1),
261  3 gvec(2*mrho1+1:), fsq, rin(1,nplane), zin(1,nplane),
262  4 imodes)
263 
264  SELECT CASE (iter)
265  CASE DEFAULT
266  gmin = min(gmin,gnorm)
267  CALL evolve (g11)
268 !
269 ! RUDIMENTARY TIME STEP CONTROL
270 !
271  IF (gnorm/gmin .gt. too_large) irst = 2
272  IF (irst.eq.2 .or. gmin.eq.gnorm) CALL restart (irst)
273 
274  IF (gnorm.lt.gtrig .and. imodes.lt.mrho) THEN
275  imodes = imodes + 1
276  CALL restart (irst)
277  irst = 2
278  delt = delt/.95
279  ENDIF
280  IF (mod(iter,nstep).ne.0 .and. gnorm.gt.ftol**2) cycle
281  CASE (2)
282  g11 = gnorm
283  cycle
284  CASE (1)
285  gmin = gnorm
286  imodes = min(4,mrho)
287  END SELECT
288 
289  gout = sqrt(gnorm)
290  modeno = imodes
291  IF (iter .eq. 1) modeno = mpol - 1
292  IF (qexp .gt. zero) specw = specw**(one/qexp)
293 CEAL PRINT 110, iter, fsq, gout, specw, MODeno, delt
294 CEAL WRITE (3, 110) iter, fsq, gout, specw, MODeno, delt
295  IF (gnorm.lt.ftol**2 .and. imodes.eq.mrho
296  1 .or. fsq.lt.ftol) EXIT
297  END DO
298 
299  110 FORMAT(i8,1p,2e16.3,0p,f10.2,i8,1p,e10.2)
300 
301 !
302 ! STORE FINAL ANSWER FOR FINAL PHI-TRANSFORM OF RHOC, RHOS
303 !
304  result(1:2*mpol,nplane) = xvec(1:2*mpol)
305 
306  END DO
307 !
308 ! OUTPUT LOOP
309 !
310 CEAL PRINT 330, pexp, qexp, mexp
311 CEAL WRITE (3, 330) pexp, qexp, mexp
312  330 FORMAT(/' ANGLE CONSTRAINTS WERE APPLIED ',/,
313  1 ' BASED ON RM**2 + ZM**2 SPECTRUM WITH P = ',f8.2,' AND Q = ',
314  2 f8.2,/,' POLAR DAMPING EXPONENT = ',i2/,' TIME: ',1pe10.2,
315  3 ' SEC.'/)
316 !
317 ! PERFORM PHI-FOURIER TRANSFORM
318 !
319  CALL fftrans (result(1,1:nphi), result(1+mpol,1:nphi),
320  1 result(2:mpol,1:nphi), result(mpol+2:2*mpol,1:nphi), rbc,
321  1 zbs, rbs, zbc, rmnaxis, zmnaxis)
322 
323  DEALLOCATE (xvec, gvec, xdot, xstore, angle)
324 
325  END SUBROUTINE scrunch
326 
327  SUBROUTINE fixaray(pexp, qexp, mexp)
328 C-----------------------------------------------
329 C M o d u l e s
330 C-----------------------------------------------
331  USE vname0
332  USE vname1
333  IMPLICIT NONE
334 C-----------------------------------------------
335 C D u m m y A r g u m e n t s
336 C-----------------------------------------------
337  INTEGER :: mexp
338  REAL(rprec) :: pexp, qexp
339 C-----------------------------------------------
340 C L o c a l V a r i a b l e s
341 C-----------------------------------------------
342  INTEGER :: l, ntor, nn0, n, m
343  REAL(rprec) :: sgn
344 C-----------------------------------------------
345 !
346 ! This routine stores toroidal and poloidal MODe number arrays
347 !
348 ! MPOL = NUMBER OF POLOIDAL MODES USED IN CURVE-FIT
349 ! NPHI = NUMBER OF EQUALLY-SPACED PHI PLANES PER FIELD
350 ! PERIOD IN WHICH THE CURVE-FITTING IS PERFORMED
351 ! IMPORTANT: NPHI MUST BE EVEN FOR NON-SYMMETRIC SYSTEMS
352 !
353 ! MPNT = NUMBER OF R AND Z MODES APPEARING IN FIT
354 ! NTHETA = NUMBER OF THETA POINTS IN EACH PHI PLANE
355 ! N2 = TOTAL NUMBER OF RESIDUALS IN FSQ (PER PHI PLANE)
356 ! == 2 (m=0 MODes) + 2*mrho (rhoc, rhos) + ntheta
357 ! NFP = NUMBER OF TOROIDAL FIELD PERIODS (IRRELEVANT)
358 !
359  mrho = mpol-1
360  nphi2 = nphi/2
361  IF (nphi2 .eq. 0) nphi2 = 1
362  n2 = 2*(mrho+1) + ntheta
363 
364  l = 0
365  ntor = max0(1,nphi - 1)
366  IF (nphi .eq. 2) ntor = 2
367 
368  nn0 = 1 - (ntor + 1)/2
369  DO n = 1, ntor
370  nn(n) = (nn0 + (n - 1))*nfp
371  END DO
372  DO m = 1, mpol
373  mm(m) = m - 1
374  DO n = 1, ntor
375  IF (mm(m).ne.0 .or. nn(n).ge.0) THEN
376  l = l + 1
377  m1(l) = mm(m)
378  n1(l) = nn(n)
379  ENDIF
380  END DO
381  END DO
382  mpnt = l
383  dnorm = 2._dbl/ntheta
384 
385  DO m = 0, mpol
386  dm1(m) = m
387  IF (m .eq. 0) THEN
388  xmpq(m,1) = zero
389  xmpq(m,2) = zero
390  ELSE
391  xmpq(m,1) = dm1(m)**(pexp+qexp)
392  xmpq(m,2) = dm1(m)**(pexp)
393  ENDIF
394  END DO
395 
396 !
397 ! FOR SECOND DERIVATIVE-TYPE CONSTRAINT, CHOOSE SGN = -one
398 !
399  sgn = one !First derivative-TYPE constraint
400 ! sgn = -one !Second derivative-TYPE constraint
401  t1m(1) = one
402  DO m = 2, mrho
403  t1m(m) = (real(m-1, rprec)/m)**mexp
404  END DO
405  DO m = 1,mrho-2
406  t2m(m) = sgn*(real(m+1, rprec)/m)**mexp
407  END DO
408  t2m(mrho-1) = zero
409  t2m(mrho) = zero
410 
411  END SUBROUTINE fixaray
412  SUBROUTINE inguess(rin, zin, angle)
413 C-----------------------------------------------
414 C M o d u l e s
415 C-----------------------------------------------
416  USE vname0
417  USE vname1
418  IMPLICIT NONE
419 C-----------------------------------------------
420 C D u m m y A r g u m e n t s
421 C-----------------------------------------------
422  REAL(rprec), DIMENSION(ntheta,nphi) :: rin, zin, angle
423 C-----------------------------------------------
424 C L o c a l V a r i a b l e s
425 C-----------------------------------------------
426  INTEGER :: i, inside, jskip, j1, j2
427 C-----------------------------------------------
428 c**********************************************************************
429 c This SUBROUTINE obtains initial guesses for the centroid at each
430 c toroidal plane. By DEFAULT, the polar axis is set equal to this
431 c centroid. It is imperative that the polar axis be enclosed by
432 c the surface (otherwise the computed theta angles will not span
433 c [0,2pi]). For certain complex cross-sections, this simple estimate
434 c of the polar axis may fail, and the USEr must supply values for
435 c raxis(i),zaxis(i). In addition, for non-starlike DOmains, the
436 c points along the curve must be monitonically increasing as a
437 c FUNCTION of the arclength from ANY fixed poINT on the curve. This
438 c ordering is attempted by the SUBROUTINE ORDER, but may fail in
439 c general IF too few theta points are USEd.
440 c
441 c Modified 8/13/97, J. Breslau. The SUBROUTINE now tries several guesses
442 c for the magnetic axis position IF the first fails. The mean position
443 c of pairs of points is taken until one is found to be inside.
444 c**********************************************************************
445 c COMPUTE CENTROID AND
446 c ORDER POINTS ON FLUX SURFACE AT EACH TOROIDAL ANGLE
447 c**********************************************************************
448 c PRINT *, 'ORDERING SURFACE POINTS'
449  r0n = 0 ; z0n = 0
450  DO i = 1, nphi
451  r0n(i) = r0n(i) + sum(rin(:,i))/real(ntheta,rprec)
452  z0n(i) = z0n(i) + sum(zin(:,i))/real(ntheta,rprec)
453  raxis(i) = r0n(i)
454  zaxis(i) = z0n(i)
455  CALL order (rin(1,i), zin(1,i), raxis(i), zaxis(i), inside)
456  IF (inside .eq. 0) THEN
457  jskip = ntheta/2
458  40 CONTINUE
459  jskip = jskip - 1
460  IF (jskip .le. 0) THEN
461  print *, 'Could not find INTernal point'
462  WRITE (3, *) 'Could not find INTernal point'
463  stop
464  ENDIF
465  j1 = 1
466  j2 = j1 + jskip
467  50 CONTINUE
468  raxis(i) = 0.5_dbl*(rin(j1,i)+rin(j2,i))
469  zaxis(i) = 0.5_dbl*(zin(j1,i)+zin(j2,i))
470  CALL order (rin(1,i), zin(1,i), raxis(i), zaxis(i), inside)
471  IF (inside .eq. 1) THEN
472  r0n(i) = raxis(i)
473  z0n(i) = zaxis(i)
474  ELSE
475  j1 = j1 + 1
476  IF (j1 .gt. ntheta) go to 40
477  j2 = j2 + 1
478  IF (j2 .gt. ntheta) j2 = 1
479  go to 50
480  ENDIF
481  ENDIF
482  END DO
483 c**********************************************************************
484 c COMPUTE OPTIMUM ANGLE OFFSETS FOR M = 1 MODES
485 c**********************************************************************
486  CALL getangle (rin, zin, angle, r0n, z0n)
487 
488  END SUBROUTINE inguess
489  SUBROUTINE amplitud(rcenter, zcenter, angin, r0c, z0c, rhoc, rhos,
490  1 xpts, xin, yin)
491 C-----------------------------------------------
492 C M o d u l e s
493 C-----------------------------------------------
494  USE vname0
495  USE vname1
496  IMPLICIT NONE
497 C-----------------------------------------------
498 C D u m m y A r g u m e n t s
499 C-----------------------------------------------
500  REAL(rprec) rcenter, zcenter, r0c, z0c
501  REAL(rprec), DIMENSION(ntheta) :: angin, xpts, xin, yin
502  REAL(rprec), DIMENSION(0:mrho-1) :: rhoc, rhos
503 C-----------------------------------------------
504 C L o c a l V a r i a b l e s
505 C-----------------------------------------------
506  INTEGER :: mrz
507  INTEGER :: m, j
508  REAL(rprec) :: xmult, arg, xi, yi, t1, t2,
509  1 r1c(mu), r1s(mu), z1c(mu), z1s(mu), tnorm
510 C-----------------------------------------------
511 c*****************************************************************
512 c This SUBROUTINE assigns initial guesses for angles and
513 c Fourier MODe amplitudes to the appropriate components of
514 c the xvec array
515 c*****************************************************************
516  r0c = rcenter
517  z0c = zcenter
518  xpts(:ntheta) = angin(:ntheta)
519 
520  mrz = mpol-1
521  xmult = 2._dbl/ntheta
522 
523  r1c = zero
524  r1s = zero
525  z1c = zero
526  z1s = zero
527  rhoc = zero
528  rhos = zero
529 
530  DO m = 1, mrz
531  DO j = 1, ntheta
532  arg = angin(j)
533  xi = xmult*(xin(j)-rcenter)
534  yi = xmult*(yin(j)-zcenter)
535  r1c(m) = r1c(m) + cos(m*arg)*xi
536  r1s(m) = r1s(m) + sin(m*arg)*xi
537  z1c(m) = z1c(m) + cos(m*arg)*yi
538  z1s(m) = z1s(m) + sin(m*arg)*yi
539  END DO
540  END DO
541 
542  r10 = sqrt( r1c(1)**2 + r1s(1)**2 + z1c(1)**2 + z1s(1)**2 )
543 CEAL WRITE(3,'(/,3(a,1pe10.3))')
544 CEAL 1 ' ]RAXIS = ', rcenter,' ZAXIS = ', zcenter,' R10 = ',r10
545 CEAL WRITE(*,'(/,3(a,1pe10.3))')
546 CEAL 1 ' RAXIS = ', rcenter,' ZAXIS = ', zcenter,' R10 = ',r10
547 
548 *
549 * INITIALIZE POLAR RHOs for m = 0 thru mrz-1
550 *
551  DO m = 0, mrz-1
552  IF( m.le.1 )then
553  t1 = t1m(m+1)
554  rhoc(m) = 0.5_dbl*(r1c(m+1) + z1s(m+1))/t1
555  rhos(m) = 0.5_dbl*(r1s(m+1) - z1c(m+1))/t1
556  ELSE
557  t1 = t1m(m+1)
558  t2 = t2m(m-1)
559  tnorm = 0.5_dbl/(t1**2 + t2**2)
560  rhoc(m) = tnorm*( (r1c(m+1) + z1s(m+1))*t1
561  1 + (r1c(m-1) - z1s(m-1))*t2 )
562  rhos(m) = tnorm*( (r1s(m+1) - z1c(m+1))*t1
563  1 + (r1s(m-1) + z1c(m-1))*t2 )
564  END IF
565  END DO
566 
567  END SUBROUTINE amplitud
568  SUBROUTINE funct(r0c, z0c, rhoc, rhos, xpts, gr0c, gz0c, grhoc,
569  1 grhos, gpts, fsq, xin, yin, mrho_in)
570 C-----------------------------------------------
571 C M o d u l e s
572 C-----------------------------------------------
573  USE vname0
574  USE vname1
575  IMPLICIT NONE
576 C-----------------------------------------------
577 C D u m m y A r g u m e n t s
578 C-----------------------------------------------
579  INTEGER mrho_in
580  REAL(rprec) fsq, r0c, z0c, gr0c, gz0c
581  REAL(rprec), DIMENSION(ntheta) :: xpts, gpts, xin, yin
582  REAL(rprec), DIMENSION(0:mrho_in-1) ::
583  1 rhoc, rhos, grhoc, grhos
584 C-----------------------------------------------
585 C L o c a l V a r i a b l e s
586 C-----------------------------------------------
587  INTEGER :: m, mrho1
588  REAL(rprec),DIMENSION(ntheta) ::
589  1 gcon, gtt, r1, z1, rt1, zt1
590  REAL(rprec), DIMENSION(ntheta,0:mrho_in) :: COSa, SINa
591  REAL(rprec) :: denom, rmc_p, zms_p, rms_p, zmc_p,
592  1 t1, t2, tnorm
593 C-----------------------------------------------
594 !
595 !
596 ! FORCES: dW/dRmn = SUM(i)[ R(i) - RIN(i)]*COS(m*u[i]) ....
597 ! dW/dZmn = SUM(i)[ Z(i) - ZIN(i)]*SIN(m*u[i]) ....
598 ! dW/du[i] = rt(i)*[R(i)-RIN(i)] + zt(i)*[Z(i) - ZIN(i)]
599 ! THE NORM ON THE ANGLE FORCE (dW/du[i]) FOLLOWS FROM NEWTONS
600 ! LAW AND IS APPROXIMATELY GTT = RT**2 + ZT**2
601 !
602  denom = zero
603  specw = zero
604 
605  mrho1 = mrho_in-1
606  r1(:ntheta) = -xin(:ntheta)
607  z1(:ntheta) = -yin(:ntheta)
608  rt1(:ntheta) = zero
609  zt1(:ntheta) = zero
610  cosa(:,0) = one
611  cosa(:,1) = cos(xpts)
612  sina(:,0) = zero
613  sina(:,1) = sin(xpts)
614 
615 !
616 ! COMPUTE CURVE R1 = R - Rin, Z1 = Z - Zin FOR INPUT POINTS
617 ! NOTE DIMENSIONS: Rm(0:MRHO), Zm(0:MRHO), but RHO(0:MRHO-1)
618 !
619  DO m = 2, mrho_in
620  cosa(:,m) = cosa(:,m-1)*cosa(:,1) - sina(:,m-1)*sina(:,1)
621  sina(:,m) = sina(:,m-1)*cosa(:,1) + cosa(:,m-1)*sina(:,1)
622  END DO
623 
624  DO m = 0, mrho_in
625  CALL getrz(rmc_p,rms_p,zmc_p,zms_p,r0c,z0c,rhoc,rhos,
626  1 m,mrho_in)
627 !
628 ! COMPUTE SPECTRAL WIDTH OF CURVE
629 !
630  t2 = rmc_p*rmc_p + zmc_p*zmc_p + rms_p*rms_p + zms_p*zms_p
631  denom = denom + t2*xmpq(m,2)
632  specw = specw + xmpq(m,1)*t2
633 
634  gtt(:ntheta) = rmc_p*cosa(:ntheta,m) + rms_p*sina(:ntheta,m)
635  gcon(:ntheta) = zmc_p*cosa(:ntheta,m) + zms_p*sina(:ntheta,m)
636  r1(:ntheta) = r1(:ntheta) + gtt(:ntheta)
637  z1(:ntheta) = z1(:ntheta) + gcon(:ntheta)
638  rt1(:ntheta) = rt1(:ntheta) + dm1(m)*(rms_p*cosa(:ntheta,m)-
639  1 rmc_p*sina(:ntheta,m))
640  zt1(:ntheta) = zt1(:ntheta) + dm1(m)*(zms_p*cosa(:ntheta,m)-
641  1 zmc_p*sina(:ntheta,m))
642  END DO
643 
644  specw = specw/denom
645 
646  gtt(:ntheta) = rt1(:ntheta)**2 + zt1(:ntheta)**2
647  gpts(:ntheta) = r1(:ntheta)*rt1(:ntheta)+z1(:ntheta)*zt1(:ntheta)
648 
649 !
650 ! COMPUTE MEAN-SQUARE DEVIATION BETWEEN POINTS AND FIT
651 !
652 ! t1 = MAXVAL(gtt(:ntheta))
653 ! gpts(:ntheta) = gpts(:ntheta)/t1
654  gpts(:ntheta) = 0.5_dbl*gpts(:ntheta)/gtt(:ntheta)
655  t1 = maxval(abs(gpts(:ntheta)))
656  t2 = 1.e-3_dbl !Arbitrary threshold on angle motion
657  IF (t1.gt.t2) THEN
658 ! Careful not to make angle points move too rapidly so they DO not cross
659  t1 = t2/t1
660  gpts(:ntheta) = t1 * gpts(:ntheta)
661  END IF
662  fsq = 0.5_dbl*dnorm*sum(r1(:ntheta)**2 + z1(:ntheta)**2)
663  fsq = sqrt(fsq)/r10
664 
665 !
666 ! COMPUTE CURVE FORCES
667 ! 1. Magnetic axis forces
668 !
669  m = 0
670  gr0c = dnorm*sum(cosa(:ntheta,m)*r1(:ntheta))
671  gz0c = dnorm*sum(cosa(:ntheta,m)*z1(:ntheta))
672 
673 !
674 ! RHO FORCES
675 ! (NOTE: RHO for m=0 IS EQUIVALENT TO THE L(v)-TERM => ARCLENGTH, Eq.14 H-B paper)
676 !
677  DO m = 0, mrho1
678  IF (m .le. 1) THEN
679  t1 = dnorm/t1m(m+1)
680  grhoc(m) = t1*sum(cosa(:ntheta,m+1)*r1(:ntheta) +
681  1 sina(:ntheta,m+1)*z1(:ntheta))
682  grhos(m) = t1*sum(sina(:ntheta,m+1)*r1(:ntheta) -
683  1 cosa(:ntheta,m+1)*z1(:ntheta))
684  ELSE
685  t1 = t1m(m+1)
686  t2 = t2m(m-1)
687  tnorm = dnorm/(t1*t1 + t2*t2)
688  t1 = t1*tnorm
689  t2 = t2*tnorm
690  grhoc(m) = sum((cosa(:ntheta,m+1)*r1(:ntheta) +
691  1 sina(:ntheta,m+1)*z1(:ntheta))*t1 +
692  2 (cosa(:ntheta,m-1)*r1(:ntheta) -
693  3 sina(:ntheta,m-1)*z1(:ntheta))*t2)
694  grhos(m) = sum((sina(:ntheta,m+1)*r1(:ntheta) -
695  1 cosa(:ntheta,m+1)*z1(:ntheta))*t1 +
696  2 (sina(:ntheta,m-1)*r1(:ntheta) +
697  3 cosa(:ntheta,m-1)*z1(:ntheta))*t2)
698  ENDIF
699  END DO
700 
701  grhos(0) = zero !Enforce constraINT on toroidal angle
702 
703  gnorm = sum(grhoc(0:mrho1)*grhoc(0:mrho1) +
704  1 grhos(0:mrho1)*grhos(0:mrho1)) +
705  2 gr0c**2 + gz0c**2
706  gnorm = gnorm/r10**2
707 
708  gnorm = gnorm + dnorm*sum(gpts(:ntheta)*gpts(:ntheta))
709 
710  END SUBROUTINE funct
711  SUBROUTINE evolve(g11)
712 C-----------------------------------------------
713 C M o d u l e s
714 C-----------------------------------------------
715  USE vname0
716  USE vname1
717  IMPLICIT NONE
718 C-----------------------------------------------
719 C D u m m y A r g u m e n t s
720 C-----------------------------------------------
721  REAL(rprec) g11
722 C-----------------------------------------------
723 C L o c a l V a r i a b l e s
724 C-----------------------------------------------
725  REAL(rprec), PARAMETER :: bmax = 0.15_dbl
726  REAL(rprec) :: ftest, dtau, otav, b1, fac
727 C-----------------------------------------------
728 
729  ftest = gnorm/g11
730  dtau = abs(one - ftest)
731  g11 = gnorm
732  otav = dtau/delt
733  dtau = delt*otav + 1.e-3_dbl
734  dtau = min(bmax,dtau)
735  b1 = one - 0.5_dbl*dtau
736  fac = one/(one + 0.5_dbl*dtau)
737  xdot(:n2) = fac*(xdot(:n2)*b1-delt*gvec(:n2))
738  xvec(:n2) = xvec(:n2) + xdot(:n2)*delt
739 
740  END SUBROUTINE evolve
741  SUBROUTINE fftrans(r0c, z0c, rhoc, rhos, rbc, zbs, rbs, zbc,
742  1 rmnaxis, zmnaxis)
743 C-----------------------------------------------
744 C M o d u l e s
745 C-----------------------------------------------
746  USE vname0
747  USE vname1
748  IMPLICIT NONE
749 C-----------------------------------------------
750 C D u m m y A r g u m e n t s
751 C-----------------------------------------------
752  REAL(rprec), DIMENSION(nphi) :: r0c, z0c
753  REAL(rprec), DIMENSION(0:mrho-1,nphi) :: rhoc, rhos
754  REAL(rprec), DIMENSION(0:mpol-1,-nphi2:nphi2) ::
755  1 rbc, zbs, rbs, zbc
756  REAL(rprec), DIMENSION(0:nphi2) :: rmnaxis, zmnaxis
757 C-----------------------------------------------
758 C L o c a l V a r i a b l e s
759 C-----------------------------------------------
760  INTEGER :: i, mn, mREAL, nreal
761  REAL(rprec), DIMENSION(nv) :: INTgrate, argi
762  REAL(rprec) :: delphi,dn,rmc_p,zms_p,rms_p,zmc_p,
763  1 arg,tcosn,tsinn
764 C-----------------------------------------------
765 !
766 ! PERFORM FOURIER TRANSFORM IN phi
767 !
768  delphi = one/nphi
769  DO i = 1, nphi
770  intgrate(i) = delphi
771  argi(i) = twopi*(i - 1)/real(nphi*nfp,rprec)
772  END DO
773 
774  rbc = 0; zbs = 0; rbs = 0; zbc = 0
775 
776  DO mn = 1, mpnt
777  mreal = m1(mn)
778  nreal = n1(mn)/nfp
779  dn = real(n1(mn))
780  DO i = 1, nphi
781  CALL getrz(rmc_p,rms_p,zmc_p,zms_p,r0c(i),z0c(i),
782  1 rhoc(0,i),rhos(0,i),mreal,mrho)
783  arg = dn*argi(i)
784  tcosn = cos(arg)
785  tsinn = sin(arg)
786  rbc(mreal,nreal) = rbc(mreal,nreal) + intgrate(i)*(tcosn*
787  1 rmc_p + tsinn*rms_p)
788  zbs(mreal,nreal) = zbs(mreal,nreal) + intgrate(i)*(tcosn*
789  1 zms_p - tsinn*zmc_p)
790  zbc(mreal,nreal) = zbc(mreal,nreal) + intgrate(i)*(tcosn*
791  1 zmc_p + tsinn*zms_p)
792  rbs(mreal,nreal) = rbs(mreal,nreal) + intgrate(i)*(tcosn*
793  1 rms_p - tsinn*rmc_p)
794  END DO
795  IF (mreal.eq.0 .and. nreal.eq.0) THEN
796  rmnaxis(0) = dot_product(intgrate(:nphi),raxis(:nphi))
797  zmnaxis(0) = zero
798  ELSE IF (mreal.eq.0 .and. nreal.gt.0) THEN
799  rbc(0,nreal) = 2*rbc(0,nreal)
800  rbs(0,nreal) = 2*rbs(0,nreal)
801  zbc(0,nreal) = 2*zbc(0,nreal)
802  zbs(0,nreal) = 2*zbs(0,nreal)
803  rmnaxis(nreal) = zero
804  zmnaxis(nreal) = zero
805  rmnaxis(nreal) = rmnaxis(nreal) + sum(2*intgrate(:nphi)*
806  1 raxis(:nphi)*cos(dn*argi(:nphi)))
807  zmnaxis(nreal) = zmnaxis(nreal) - sum(2*intgrate(:nphi)*
808  1 zaxis(:nphi)*sin(dn*argi(:nphi)))
809  ENDIF
810  END DO
811 
812  END SUBROUTINE fftrans
813  SUBROUTINE printit(rin,zin,rbc,zbs,rbs,zbc,rmnaxis,zmnaxis,g2)
814 C-----------------------------------------------
815 C M o d u l e s
816 C-----------------------------------------------
817  USE vname0
818  USE vname1
819  USE mapout
820  IMPLICIT NONE
821 C-----------------------------------------------
822 C D u m m y A r g u m e n t s
823 C-----------------------------------------------
824  REAL(rprec), DIMENSION(*) :: rin, zin
825  REAL(rprec), DIMENSION(0:mpol - 1,-nphi2:nphi2) ::
826  1 rbc, zbs, rbs, zbc
827  REAL(rprec), DIMENSION(0:nphi2) :: rmnaxis, zmnaxis
828 C-----------------------------------------------
829 C L o c a l V a r i a b l e s
830 C-----------------------------------------------
831  INTEGER :: i, m, n, n11
832  REAL(rprec) :: tol
833  CHARACTER*250 :: form_string
834  TYPE(RSZS) :: g2
835 C-----------------------------------------------
836  g2%rbc=0 ; g2%rbs=0 ; g2%zbs=0 ; g2%zbc=0 ;
837  g2%rbc(0:mu-1,0)=rbc(0:mu-1,0)
838  g2%rbs(0:mu-1,0)=rbs(0:mu-1,0)
839  g2%zbc(0:mu-1,0)=zbc(0:mu-1,0)
840  g2%zbs(0:mu-1,0)=zbs(0:mu-1,0)
841 
842  OPEN(unit=10, file='plotout', status='unknown')
843  WRITE (10, 1990) mpol, ntheta, nphi, mpol - 1, nphi2, nfp, mpnt
844  1990 FORMAT(7i10)
845  DO i = 1, nphi*ntheta
846  WRITE (10, 1995) rin(i), zin(i)
847  END DO
848  1995 FORMAT(1p,2e12.4)
849 
850 c**********************************************************************
851 c This SUBROUTINE PRINTs out data INTo the file "outcurve"
852 c The FORMAT of the MODes is compatible with input to VMEC
853 c**********************************************************************
854 CEAL WRITE (3, 10)
855  10 FORMAT(/' MB NB RBC RBS ',
856  1 'ZBC ZBS RAXIS ZAXIS')
857  tol = 1.e-6_dbl*abs(rbc(1,0))
858  DO m = 0, mpol - 1
859  DO n = -nphi2, nphi2
860 CEAL WRITE (10, 2000) rbc(m,n), zbs(m,n), rbs(m,n), zbc(m,n)
861  IF (.not.(abs(rbc(m,n)).lt.tol .and. abs(zbs(m,n)).lt.tol
862  1 .and. abs(rbs(m,n)).lt.tol .and. abs(zbc(m,n)).lt.tol))
863  2 THEN
864  IF (m.eq.0 .and. n.ge.0) THEN
865 CEAL WRITE (3, 30) m, n, rbc(m,n), rbs(m,n), zbc(m,n),
866 CEAL 1 zbs(m,n), rmnaxis(n), zmnaxis(n)
867  ELSE
868 CEAL WRITE (3, 40) m, n, rbc(m,n), rbs(m,n), zbc(m,n),
869 CEAL 1 zbs(m,n)
870  ENDIF
871  ENDIF
872  END DO
873  END DO
874  30 FORMAT(i5,i4,1p,6e12.4)
875  40 FORMAT(i5,i4,1p,4e12.4)
876  2000 FORMAT(1p,4e12.4)
877 
878 c**********************************************************************
879 c WRITE OUT IN FORMAT THAT CAN BE EXTRACTED INTO VMEC
880 c**********************************************************************
881 CEAL WRITE(3, *)
882  WRITE(67,fmt='('' MPOL = '',i2.2)')mpol
883  DO m = 0, mpol-1
884  DO n = -nphi2, nphi2
885  IF (.not.(abs(rbc(m,n)).lt.tol .and. abs(zbs(m,n)).lt.tol
886  1 .and. abs(rbs(m,n)).lt.tol .and. abs(zbc(m,n)).lt.tol))
887  2 THEN
888  n11 = nphi2/10 + 1
889  IF (n .lt. 0) n11 = n11+1
890  WRITE(form_string,'(a,4(a,i1,a,i1,a))')"(2x,",
891  1"'RBC(',i",n11,",',',i",m/10+1,",') = ',1p,e14.6,3x,",
892  2"'RBS(',i",n11,",',',i",m/10+1,",') = ',e14.6,3x,",
893  3"'ZBC(',i",n11,",',',i",m/10+1,",') = ',e14.6,3x,",
894  4"'ZBS(',i",n11,",',',i",m/10+1,",') = ',e14.6)"
895  WRITE(form_string,'(a,4(a,i1,a,i1,a))')"(2x,",
896  1"'RBC(',i",n11,",',',i",m/10+1,",') = ',1p,e14.6,3x,",
897  2"'RBS(',i",n11,",',',i",m/10+1,",') = ',e14.6,3x,/,4x,",
898  3"'ZBC(',i",n11,",',',i",m/10+1,",') = ',e14.6,3x,",
899  4"'ZBS(',i",n11,",',',i",m/10+1,",') = ',e14.6)"
900 
901  WRITE(67, form_string) n,m,rbc(m,n), n,m,rbs(m,n),
902  1 n,m,zbc(m,n), n,m,zbs(m,n)
903 c IF (m.eq.0 .and. n.ge.0) THEN
904 c WRITE (3, 30) m, n, rbc(m,n), rbs(m,n), zbc(m,n),
905 c 1 zbs(m,n), rmnaxis(n), zmnaxis(n)
906  END IF
907  END DO
908  END DO
909 
910  END SUBROUTINE printit
911  SUBROUTINE getrz(rmc,rms,zmc,zms,r0c,z0c,rhoc,rhos,m,mrho_in)
912  USE vname1
913 C-----------------------------------------------
914 C D u m m y A r g u m e n t s
915 C-----------------------------------------------
916  INTEGER m, mrho_in
917  REAL(rprec), DIMENSION(0:mrho_in-1) :: rhoc, rhos
918  REAL(rprec) :: rmc, rms, zmc, zms, r0c, z0c
919 C-----------------------------------------------
920 C L o c a l V a r i a b l e s
921 C-----------------------------------------------
922  INTEGER :: mrho1
923 C-----------------------------------------------
924 
925  mrho1 = mrho_in - 1
926  rhos(0) = zero !Enforce ConstraINT on toroidal angle
927 
928  IF (m .eq. 0) THEN
929  rmc = r0c
930  zmc = z0c
931  rms = zero
932  zms = zero
933  ELSE IF (m .lt. mrho1) THEN
934  rmc = (t1m(m)*rhoc(m-1) + t2m(m)*rhoc(m+1))
935  zms = (t1m(m)*rhoc(m-1) - t2m(m)*rhoc(m+1))
936  rms = (t1m(m)*rhos(m-1) + t2m(m)*rhos(m+1))
937  zmc =-(t1m(m)*rhos(m-1) - t2m(m)*rhos(m+1))
938  ELSE !Can change highest m constraints here...
939  rmc = t1m(m)*rhoc(m-1) * hb_parameter
940  zms = rmc * hb_parameter
941  rms = t1m(m)*rhos(m-1) * hb_parameter
942  zmc =-rms * hb_parameter
943  ENDIF
944 
945  END SUBROUTINE getrz
946  SUBROUTINE getangle(rval, zval, angle, rcenter, zcenter)
947 C-----------------------------------------------
948 C M o d u l e s
949 C-----------------------------------------------
950  USE vname0
951  USE vname1
952  IMPLICIT NONE
953 C-----------------------------------------------
954 C D u m m y A r g u m e n t s
955 C-----------------------------------------------
956  REAL(rprec), DIMENSION(ntheta, nphi) :: rval, zval, angle
957  REAL(rprec), DIMENSION(nv) :: rcenter, zcenter
958 C-----------------------------------------------
959 C L o c a l V a r i a b l e s
960 C-----------------------------------------------
961  INTEGER :: i, j, iterate
962  REAL(rprec), DIMENSION(nv) ::
963  1 rcos, rsin, zcos, zsin, phiangle
964  REAL(rprec) :: xc, yc, dnum, denom, delangle
965 C-----------------------------------------------
966 c**********************************************************************
967 c Compute angle offset consistent with constraINT Z1n = Z1,-n
968 c Note: This is DOne iteratively, SINce elongation is unknown
969 c**********************************************************************
970  DO i = 1, nphi
971  DO j = 1, ntheta
972  angle(j,i) = twopi*(j - 1)/real(ntheta,rprec)
973  END DO
974  END DO
975  DO iterate = 1, 5
976  DO i = 1, nphi
977  rcos(i) = zero
978  rsin(i) = zero
979  zcos(i) = zero
980  zsin(i) = zero
981  DO j = 1, ntheta
982  xc = rval(j,i) - rcenter(i)
983  yc = zval(j,i) - zcenter(i)
984  rcos(i) = rcos(i) + cos(angle(j,i))*xc
985  rsin(i) = rsin(i) + sin(angle(j,i))*xc
986  zcos(i) = zcos(i) + cos(angle(j,i))*yc
987  zsin(i) = zsin(i) + sin(angle(j,i))*yc
988  END DO
989  END DO
990 c**********************************************************************
991 c Compute new angles starting from offset phiangle(i)
992 c**********************************************************************
993  dnum = zero
994  denom = zero
995  dnum = sum(zsin(:nphi))
996  denom = sum(rcos(:nphi))
997  IF (denom .ne. zero) THEN
998  elongate = dnum/denom
999  ELSE
1000  elongate = 1.e10
1001  END IF
1002 
1003  delangle = zero
1004  DO i = 1, nphi
1005  phiangle(i) = atan2(elongate*zcos(i)-rsin(i),
1006  1 elongate*zsin(i)+rcos(i))
1007  delangle = max(delangle,abs(phiangle(i)))
1008  angle(:,i) = angle(:,i) + phiangle(i)
1009  END DO
1010  IF (delangle < 0.02_dbl) EXIT
1011  END DO
1012 CEAL WRITE (*, 1010) elongate, raxis(1), zaxis(1)
1013 CEAL WRITE (3, 1010) elongate, raxis(1), zaxis(1)
1014 CEAL WRITE (*, 1020) ntheta, nphi
1015 CEAL WRITE (3, 1020) ntheta, nphi
1016  1010 FORMAT(' Average elongation = ',1pe10.4,/,' Raxis = ',1pe12.4,
1017  1 ' Zaxis = ',1pe12.4)
1018  1020 FORMAT(' Number of Theta Points Matched = ',i5,/,
1019  1 ' Number Phi Planes = ',i5)
1020 
1021  END SUBROUTINE getangle
1022  SUBROUTINE restart(irst)
1023 C-----------------------------------------------
1024 C M o d u l e s
1025 C-----------------------------------------------
1026  USE vname0
1027  USE vname1
1028  IMPLICIT NONE
1029 C-----------------------------------------------
1030 C D u m m y A r g u m e n t s
1031 C-----------------------------------------------
1032  INTEGER irst
1033 C-----------------------------------------------
1034 C L o c a l V a r i a b l e s
1035 C-----------------------------------------------
1036  INTEGER, SAVE :: nresets = 0
1037 C-----------------------------------------------
1038 !
1039 ! This routine either stores an accepted value of the local solution
1040 ! (irst = 1) or reset the PRESENT solution to a previous value (irst = 2)
1041 !
1042  SELECT CASE (irst)
1043 
1044  CASE DEFAULT
1045  xstore(:n2) = xvec(:n2)
1046  RETURN
1047 
1048  CASE (2)
1049  xdot(:n2) = zero
1050  xvec(:n2) = xstore(:n2)
1051  delt = .95*delt
1052  irst = 1
1053  nresets = nresets + 1
1054  IF (nresets .ge. 100) THEN
1055  print *, ' Time step reduced 100 times without convergence'
1056  stop
1057  ENDIF
1058  RETURN
1059  END SELECT
1060 
1061  END SUBROUTINE restart
1062  SUBROUTINE order(rval, zval, xaxis, yaxis, inside)
1063 C-----------------------------------------------
1064 C M o d u l e s
1065 C-----------------------------------------------
1066  USE vname0
1067  USE vname1
1068  IMPLICIT NONE
1069 C-----------------------------------------------
1070 C D u m m y A r g u m e n t s
1071 C-----------------------------------------------
1072  INTEGER inside
1073  REAL(rprec) xaxis, yaxis
1074  REAL(rprec), DIMENSION(*) :: rval, zval
1075 C-----------------------------------------------
1076 C L o c a l V a r i a b l e s
1077 C-----------------------------------------------
1078  INTEGER :: i, ip1, i1, j, next
1079  REAL(rprec), DIMENSION(nu) :: tempr, tempz
1080  REAL(rprec) ::
1081  1 newdist,olddist,shortest,saver,savez,residue,x,y,dx,dy
1082 C-----------------------------------------------
1083 c**********************************************************************
1084 c Program ORDER : orders points on a magnetic surface at a
1085 c fixed toroidal plane and assigns right-handed circulation
1086 c around flux surface. XAXIS, YAXIS: Polar-TYPE axis (must lie
1087 c inside curve to check SIGN of rotation)
1088 c
1089 c Modified 8/13/97, J. Breslau. No longer quits IF the axis is not
1090 c contained by the points. Instead, RETURNs inside=0 so inguess can
1091 c try again.
1092 c**********************************************************************
1093  olddist = 1.e20_dbl
1094  DO i = 1, ntheta - 1
1095  ip1 = i + 1
1096  i1 = i
1097  shortest = 1.e20_dbl
1098  15 CONTINUE
1099  DO j = ip1, ntheta
1100  IF (i1 > 1) olddist = (rval(i1-1)-rval(j))**2 + (zval(i1-1)-
1101  1 zval(j))**2
1102  newdist = (rval(i1)-rval(j))**2 + (zval(i1)-zval(j))**2
1103  IF (newdist.le.olddist .and. newdist<shortest) THEN
1104  next = j
1105  shortest = newdist
1106  ENDIF
1107  END DO
1108 c**********************************************************************
1109 c Swap nearest poINT (next) with current poINT (ip1)
1110 c**********************************************************************
1111  IF (shortest .ge. 1.e10_dbl) THEN
1112  saver = rval(i-1)
1113  rval(i-1) = rval(i)
1114  rval(i) = saver
1115  savez = zval(i-1)
1116  zval(i-1) = zval(i)
1117  zval(i) = savez
1118  i1 = i1 - 1
1119  ip1 = ip1 - 1
1120  go to 15
1121  ENDIF
1122  saver = rval(ip1)
1123  rval(ip1) = rval(next)
1124  rval(next) = saver
1125  savez = zval(ip1)
1126  zval(ip1) = zval(next)
1127  zval(next) = savez
1128  END DO
1129 c**********************************************************************
1130 c Check that xaxis,yaxis is inside surface and
1131 c ascertain that the angle rotates counterclockwise
1132 c using Cauchy theorem in "complex"-plane
1133 c**********************************************************************
1134  inside = 1
1135  residue = zero
1136  DO i = 1, ntheta - 1
1137  x = 0.5_dbl*(rval(i)+rval(i+1)) - xaxis
1138  y = 0.5_dbl*(zval(i)+zval(i+1)) - yaxis
1139  dx = rval(i+1) - rval(i)
1140  dy = zval(i+1) - zval(i)
1141  residue = residue + (x*dy - y*dx)/(x**2 + y**2 + 1.e-10_dbl)
1142  END DO
1143  x = 0.5_dbl*(rval(1)+rval(ntheta)) - xaxis
1144  y = 0.5_dbl*(zval(1)+zval(ntheta)) - yaxis
1145  dx = rval(1) - rval(ntheta)
1146  dy = zval(1) - zval(ntheta)
1147  residue = residue + (x*dy - y*dx)/(x**2 + y**2 + 1.e-10_dbl)
1148 
1149  IF (residue < (-0.9_dbl*twopi)) THEN
1150  DO i = 2, ntheta
1151  j = ntheta - i + 2
1152  tempr(i) = rval(j)
1153  tempz(i) = zval(j)
1154  END DO
1155  rval(2:ntheta) = tempr(2:ntheta)
1156  zval(2:ntheta) = tempz(2:ntheta)
1157  ELSE IF (abs(residue) < 0.9_dbl*twopi) THEN
1158  print *, ' mag. axis not enclosed by bndry; trying again'
1159  WRITE (3, *) ' mag. axis not enclosed by bndry; trying again'
1160  inside = 0
1161 c STOP
1162  ENDIF
1163 
1164  END SUBROUTINE order
1165  SUBROUTINE evals_boundary(g2,xp,yp,icount)
1166  USE vname0
1167  USE mapout
1168  REAL(rprec), DIMENSION(0:mu-1,0:0) :: rbc,zbs,rbs,zbc
1169  real*4, DIMENSION(:), ALLOCATABLE :: rb, zb
1170  real*4, DIMENSION(*) :: xp, yp
1171  TYPE(RSZS) :: g2
1172  mpol=mu-1
1173  ntor=0
1174  ntheta=g2%nthet
1175  rbc(0:mu-1,0)=g2%rbc(0:mu-1,0)
1176  rbs(0:mu-1,0)=g2%rbs(0:mu-1,0)
1177  zbc(0:mu-1,0)=g2%zbc(0:mu-1,0)
1178  zbs(0:mu-1,0)=g2%zbs(0:mu-1,0)
1179  pi=atan(1.d0)*4.d0
1180  pit = 2.d0*pi/(ntheta - 1)
1181  v=0
1182  IF( ALLOCATED(rb) ) DEALLOCATE(rb,zb)
1183  ALLOCATE ( rb(ntheta+1), zb(ntheta+1) )
1184  icount=0
1185  DO j=1,ntheta,2
1186  u=pit*(j-1)
1187  r=0
1188  z=0
1189  n=0
1190  DO m=0,mpol
1191 ! DO n=-ntor,ntor
1192  xxm=m
1193  arg=xxm*u
1194  cosa=cos(arg);sina=sin(arg)
1195  r=r+rbc(m,n)*cosa+rbs(m,n)*sina
1196  z=z+zbc(m,n)*cosa+zbs(m,n)*sina
1197 ! ENDdo
1198  ENDdo
1199  icount=icount+1
1200  rb(icount)=r
1201  zb(icount)=z
1202  ENDdo
1203  xp(1:icount)=rb(1:icount)
1204  yp(1:icount)=zb(1:icount)
1205  RETURN
1206  END SUBROUTINE evals_boundary
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17