V3FIT
scrunch.f
1  SUBROUTINE scrunch(rin, zin, rbc, zbs, rbs, zbc,
2  1 rmnaxis, zmnaxis, ftol, niter, nstep, mexp, mntot)
3 C-----------------------------------------------
4 C M o d u l e s
5 C-----------------------------------------------
6  USE vname0
7  USE vname1
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 :: niter, nstep, mexp, mntot
13  REAL(rprec) :: ftol
14  REAL(rprec), DIMENSION(ntheta,nphi) :: rin, zin
15  REAL(rprec), DIMENSION(mntot) :: rbc, zbs, rbs, zbc
16  REAL(rprec), DIMENSION(*) :: rmnaxis, zmnaxis
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 :: pexp = 4, qexp = 1
21 C-----------------------------------------------
22 C L o c a l V a r i a b l e s
23 C-----------------------------------------------
24  INTEGER :: irst=1, nplane, imodes, iter, modeno, mrho1
25  REAL(rprec), ALLOCATABLE, DIMENSION(:,:) :: angle
26  REAL(rprec), DIMENSION(2*mpol,nphi) :: result1
27  REAL(rprec) :: too_large,
28  1 gtrig, fsq, gmin, g11, gout, time_on, time_off
29 C-----------------------------------------------
30 c
31 c PARAMETER VALUES TO BE SET IN NAME1.INC:
32 c
33 c MU: MAXIMUM NO. OF POLOIDAL MODES IN FIT
34 c NU: MAXIMUM NO. OF THETA (POLOIDAL) POINTS ALLOWED IN FIT
35 c NV: MAXIMUM NO. OF TOROIDAL PLANES IN FIT
36 c
37 !
38 
39 !
40 ! CHECK THAT INPUT DIMENSIONS ARE CONSISTENT WITH NAME1 PARAMETERS
41 !
42  deltf = 1
43  IF (mod(nphi,2).ne.0 .and. nphi.ne.1 .or. nphi.eq.0) THEN
44  print *, ' NPHI > 0 must be EVEN for non-symmetric systems'
45  stop
46  ENDIF
47  IF (ntheta > nu) THEN
48  print *, ' NTHETA input to SCRUNCH is too large'
49  stop
50  ENDIF
51  IF (nphi > nv) THEN
52  print *, ' NPHI input to SCRUNCH is too large'
53  stop
54  ENDIF
55  IF (mpol > mu) THEN
56  print *, ' MPOL input to SCRUNCH is too large'
57  stop
58  ENDIF
59  IF (nfp .le. 0) THEN
60  print *, ' NFP > 0 must be positive and exceed zero'
61  stop
62  ENDIF
63 
64 
65 !
66 ! INITIALIZE FIXED m,n ARRAYS
67 !
68  OPEN(unit=3, file='outcurve', status='unknown')
69  twopi = 8*atan(one)
70  CALL fixaray (pexp, qexp, mexp)
71 
72 !
73 ! COMPUTE INITIAL GUESSES (MUST ORDER ANGLES MONOTONICALLY FOR
74 ! NON-STARLIKE DOMAINS)
75 !
76  ALLOCATE (angle(ntheta, nphi))
77  CALL inguess (rin, zin, angle)
78 
79  too_large = 1.e6_dp
80  gtrig = 1.e-4_dp
81  mrho1 = mrho+1
82 
83  ALLOCATE (xvec(n2), gvec(n2), xdot(n2), xstore(n2))
84 !
85 ! BEGIN MAIN INTERATION LOOP
86 !
87  CALL second0(time_on)
88  DO nplane = 1, nphi
89 !
90 ! INITIALIZE M = 0 and M = 1 MODE AMPLITUDES
91 !
92 ! STACKING OF XVEC ARRAY (FOR FIXED TOROIDAL PLANE)
93 ! XVEC(1) : R0c (Symmetric components)
94 ! XVEC(2,mrho1) : Rhoc (Symmetric components)
95 ! XVEC(1+mrho1) : Z0c (Asymmetric components)
96 ! XVEC(2+mrho1,2*mrho1): Rhos (Asymmetric components)
97 ! XVEC(2*mrho1+1,n2) : Theta angle
98 !
99  xstore(:n2) = zero
100  xdot(:n2) = zero
101  xvec(:n2) = zero
102 
103  WRITE (3, 10) nplane
104  WRITE (*, 10) nplane
105  10 FORMAT(/' Fitting toroidal plane # ',i3)
106 
107  CALL amplitud (r0n(nplane), z0n(nplane), angle(1,nplane),
108  1 xvec, xvec(1+mrho1), xvec(2:mrho1), xvec(2+mrho1:2*mrho1),
109  2 xvec(2*mrho1+1:), rin(1,nplane), zin(1,nplane))
110 
111  WRITE (3, 20)
112  WRITE (*, 20)
113  20 FORMAT(/' ITERATIONS RMS ERROR FORCE GRADIENT <M>',
114  1 ' MAX m DELT')
115 
116  imodes = mrho
117  delt = deltf
118  CALL restart(0)
119 
120  DO iter = 1, niter
121  gvec(:n2) = zero
122 
123  CALL funct (xvec, xvec(1+mrho1), xvec(2:mrho1),
124  1 xvec(2+mrho1:2*mrho1),xvec(2*mrho1+1:), gvec,
125  2 gvec(1+mrho1), gvec(2:mrho1), gvec(2+mrho1:2*mrho1),
126  3 gvec(2*mrho1+1:), fsq, rin(1,nplane), zin(1,nplane),
127  4 imodes)
128 
129  SELECT CASE (iter)
130  CASE DEFAULT
131  gmin = min(gmin,gnorm)
132  CALL evolve (g11)
133 !
134 ! RUDIMENTARY TIME STEP CONTROL
135 !
136  IF (gnorm/gmin .gt. too_large) irst = 2
137  IF (irst.eq.2 .or. gmin.eq.gnorm) THEN
138  CALL restart (irst)
139  IF (irst .eq. 4) THEN
140  gnorm=gmin
141  EXIT
142  END IF
143  END IF
144 
145  IF (gnorm.lt.gtrig .and. imodes.lt.mrho) THEN
146  imodes = imodes + 1
147  CALL restart (irst)
148  irst = 2
149  delt = min(0.98_dp, delt/.975)
150  ENDIF
151  IF (mod(iter,nstep).ne.0 .and. gnorm.gt.ftol**2) cycle
152  CASE (2)
153  g11 = gnorm
154  cycle
155  CASE (1)
156  gmin = gnorm
157  imodes = min(4,mrho)
158  END SELECT
159 
160  gout = sqrt(gnorm)
161  modeno = imodes
162  IF (iter .eq. 1) modeno = mpol - 1
163  IF (qexp .gt. zero) specw = specw**(one/qexp)
164  print 110, iter, fsq, gout, specw, modeno, delt
165  WRITE (3, 110) iter, fsq, gout, specw, modeno, delt
166  IF (gnorm.lt.ftol**2 .and. imodes.eq.mrho
167  1 .or. fsq.lt.ftol) EXIT
168  END DO
169 
170  110 FORMAT(i8,1p,2e16.3,0p,f10.2,i8,1p,e10.2)
171 
172 !
173 ! STORE FINAL ANSWER FOR FINAL PHI-TRANSFORM OF RHOC, RHOS
174 !
175  result1(1:2*mpol,nplane) = xvec(1:2*mpol)
176 
177  END DO
178  CALL second0(time_off)
179 !
180 ! OUTPUT LOOP
181 !
182  print 330, mexp, pexp, qexp, time_off-time_on
183  WRITE (3, 330) mexp, pexp, qexp, time_off-time_on
184  330 FORMAT(/' ANGLE CONSTRAINTS WERE APPLIED',
185  1 ' FOR POLAR DAMPING EXPONENT (PEXP) = ',i2,/,
186  2 ' RM**2 + ZM**2 SPECTRUM COMPUTED WITH P = ',f8.2,
187  3 ' AND Q = ',f8.2,/,' TIME: ',1p,e10.2,
188  3 ' SEC.'/)
189 !
190 ! PERFORM PHI-FOURIER TRANSFORM
191 !
192  CALL fftrans (result1(1,1:nphi), result1(1+mpol,1:nphi),
193  1 result1(2:mpol,1:nphi), result1(mpol+2:2*mpol,1:nphi), rbc,
194  1 zbs, rbs, zbc, rmnaxis, zmnaxis)
195 
196  DEALLOCATE (xvec, gvec, xdot, xstore, angle)
197 
198  END SUBROUTINE scrunch
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17