1 SUBROUTINE scrunch(rin, zin, rbc, zbs, rbs, zbc,
2 1 rmnaxis, zmnaxis, ftol, niter, nstep, mexp, mntot)
12 INTEGER :: niter, nstep, mexp, mntot
14 REAL(rprec),
DIMENSION(ntheta,nphi) :: rin, zin
15 REAL(rprec),
DIMENSION(mntot) :: rbc, zbs, rbs, zbc
16 REAL(rprec),
DIMENSION(*) :: rmnaxis, zmnaxis
20 REAL(rprec),
PARAMETER :: pexp = 4, qexp = 1
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
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'
48 print *,
' NTHETA input to SCRUNCH is too large'
52 print *,
' NPHI input to SCRUNCH is too large'
56 print *,
' MPOL input to SCRUNCH is too large'
60 print *,
' NFP > 0 must be positive and exceed zero'
68 OPEN(unit=3, file=
'outcurve', status=
'unknown')
70 CALL fixaray (pexp, qexp, mexp)
76 ALLOCATE (angle(ntheta, nphi))
77 CALL inguess (rin, zin, angle)
83 ALLOCATE (xvec(n2), gvec(n2), xdot(n2), xstore(n2))
105 10
FORMAT(/
' Fitting toroidal plane # ',i3)
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))
113 20
FORMAT(/
' ITERATIONS RMS ERROR FORCE GRADIENT <M>',
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),
131 gmin = min(gmin,gnorm)
136 IF (gnorm/gmin .gt. too_large) irst = 2
137 IF (irst.eq.2 .or. gmin.eq.gnorm)
THEN
139 IF (irst .eq. 4)
THEN
145 IF (gnorm.lt.gtrig .and. imodes.lt.mrho)
THEN
149 delt = min(0.98_dp, delt/.975)
151 IF (mod(iter,nstep).ne.0 .and. gnorm.gt.ftol**2) cycle
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
170 110
FORMAT(i8,1p,2e16.3,0p,f10.2,i8,1p,e10.2)
175 result1(1:2*mpol,nplane) = xvec(1:2*mpol)
178 CALL second0(time_off)
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,
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)
196 DEALLOCATE (xvec, gvec, xdot, xstore, angle)
198 END SUBROUTINE scrunch