3 USE stel_constants,
ONLY: twopi, one, zero
28 SUBROUTINE flx2cyl(rzl_array, c_flux, r_cyl, ns, ntor,
29 1 mpol, ntmax, lthreed, lasym, iflag,
30 2 mscale, nscale, Ru, Rv, Zu, Zv)
35 INTEGER,
INTENT(out) :: iflag
36 INTEGER,
INTENT(in) :: ns, ntmax, mpol, ntor
37 LOGICAL :: lthreed, lasym
38 REAL(rprec),
DIMENSION(ns,0:ntor,0:mpol-1,3*ntmax),
39 1
INTENT(in) :: rzl_array
40 REAL(rprec),
INTENT(in) :: c_flux(3)
41 REAL(rprec),
INTENT(out) :: r_cyl(3)
42 REAL(rprec),
INTENT(in),
OPTIONAL ::
43 1 mscale(0:mpol-1), nscale(0:ntor)
44 REAL(rprec),
INTENT(out),
OPTIONAL :: Ru, Rv, Zu, Zv
48 INTEGER,
PARAMETER :: rcc = 1
52 INTEGER :: rss, rsc, rcs, zsc, zcs, zcc, zss
53 INTEGER :: istat, jslo, jshi, mpol1, m, n
54 REAL(rprec),
DIMENSION(0:ntor,0:mpol-1) ::
55 1 rmncc, rmnss, zmncs, zmnsc,
56 2 rmncs, rmnsc, zmncc, zmnss
57 REAL(rprec) :: wlo, whi, wlo_odd, whi_odd, hs1,
58 1 si, ui, vi, r11, z11
59 REAL(rprec) :: wmins(0:ntor,0:mpol-1),
60 1 wplus(0:ntor,0:mpol-1)
61 REAL(rprec) :: cosu, sinu, cosv, sinv,
62 1 cosmu(0:mpol-1), sinmu(0:mpol-1),
63 2 cosnv(0:ntor), sinnv(0:ntor),
64 3 cosnvn(0:ntor), sinnvn(0:ntor)
65 REAL(rprec) :: work1(0:mpol-1,12)
66 LOGICAL :: lru, lrv, lzu, lzv
91 si = c_flux(1); ui = c_flux(2); vi = c_flux(3)
95 IF (si .lt. zero)
THEN
96 WRITE(6, *)
' In flx2cyl, s(flux) must be > 0'
100 lru =
PRESENT(ru); lrv =
PRESENT(rv)
101 lzu =
PRESENT(zu); lzv =
PRESENT(zv)
102 IF (lrv .and. .not. lthreed) rv = 0
103 IF (lzv .and. .not. lthreed) zv = 0
112 IF (si .le. one)
THEN
116 wlo = (hs1*(jshi-1) - si)/hs1
118 IF (jslo .eq. ns) jshi = jslo
123 whi_odd = whi*sqrt(si/(hs1*(jshi-1)))
124 IF (jslo .ne. 1)
THEN
125 wlo_odd = wlo*sqrt(si/(hs1*(jslo-1)))
128 whi_odd = sqrt(si/(hs1*(jshi-1)))
135 wlo = -(si - 1)/hs1; wlo_odd = wlo
136 whi = 1 - wlo; whi_odd = whi
143 wmins(:,0:mpol1:2) = wlo
144 wplus(:,0:mpol1:2) = whi
145 wmins(:,1:mpol1:2) = wlo_odd
146 wplus(:,1:mpol1:2) = whi_odd
152 IF (ntmax .ne. 2) stop
'ntmax != 2 in flx2cyl!'
155 IF (ntmax .ne. 1) stop
'ntmax != 1 in flx2cyl!'
159 IF (ntmax .ne. 4) stop
'ntmax != 4 in flx2cyl!'
160 rss = 2; rsc = 3; rcs = 4
161 zcs = 2; zcc = 3; zss = 4
163 IF (ntmax .ne. 2) stop
'ntmax != 2 in flx2cyl!'
168 zsc = 1+ntmax; zcs = zcs+ntmax; zcc = zcc+ntmax; zss = zss+ntmax
170 rmncc = wmins*rzl_array(jslo,:,:,rcc)
171 1 + wplus*rzl_array(jshi,:,:,rcc)
172 zmnsc = wmins*rzl_array(jslo,:,:,zsc)
173 1 + wplus*rzl_array(jshi,:,:,zsc)
176 rmnss = wmins*rzl_array(jslo,:,:,rss)
177 1 + wplus*rzl_array(jshi,:,:,rss)
178 zmncs = wmins*rzl_array(jslo,:,:,zcs)
179 1 + wplus*rzl_array(jshi,:,:,zcs)
185 cosu = cos(ui); sinu = sin(ui)
186 cosv = cos(vi); sinv = sin(vi)
188 cosmu(0) = 1; sinmu(0) = 0
189 cosnv(0) = 1; sinnv(0) = 0
191 cosmu(m) = cosmu(m-1)*cosu - sinmu(m-1)*sinu
192 sinmu(m) = sinmu(m-1)*cosu + cosmu(m-1)*sinu
195 IF (
PRESENT(mscale))
THEN
196 cosmu = cosmu*mscale; sinmu = sinmu*mscale
200 cosnv(n) = cosnv(n-1)*cosv - sinnv(n-1)*sinv
201 sinnv(n) = sinnv(n-1)*cosv + cosnv(n-1)*sinv
204 IF (
PRESENT(nscale))
THEN
205 cosnv = cosnv*nscale; sinnv = sinnv*nscale
208 IF (lrv .or. lzv)
THEN
210 cosnvn(n) = n*cosnv(n)
211 sinnvn(n) =-n*sinnv(n)
224 work1(m,1) = sum(rmncc(:,m)*cosnv(:))
225 work1(m,2) = sum(zmnsc(:,m)*cosnv(:))
226 IF (lru) work1(m,3) =-m*work1(m,1)
227 IF (lzu) work1(m,4) = m*work1(m,2)
229 IF (lrv) work1(m,5) = sum(rmncc(:,m)*sinnvn(:))
230 IF (lzv) work1(m,6) = sum(zmnsc(:,m)*sinnvn(:))
231 work1(m,7) = sum(rmnss(:,m)*sinnv(:))
232 work1(m,8) = sum(zmncs(:,m)*sinnv(:))
233 IF (lru) work1(m,9) = m*work1(m,7)
234 IF (lzu) work1(m,10) =-m*work1(m,8)
235 IF (lrv) work1(m,11) = sum(rmnss(:,m)*cosnvn(:))
236 IF (lzv) work1(m,12) = sum(zmncs(:,m)*cosnvn(:))
245 r11 = sum(work1(:,1)*cosmu(:) + work1(:,7)*sinmu(:))
246 z11 = sum(work1(:,2)*sinmu(:) + work1(:,8)*cosmu(:))
247 IF (lru) ru = sum(work1(:,3)*sinmu(:) + work1(:,9)*cosmu(:))
248 IF (lzu) zu = sum(work1(:,4)*cosmu(:) + work1(:,10)*sinmu(:))
249 IF (lrv) rv = sum(work1(:,5)*cosmu(:) + work1(:,11)*sinmu(:))
250 IF (lzv) zv = sum(work1(:,6)*sinmu(:) + work1(:,12)*cosmu(:))
252 r11 = sum(work1(:,1)*cosmu(:))
253 z11 = sum(work1(:,2)*sinmu(:))
254 IF (lru) ru = sum(work1(:,3)*sinmu(:))
255 IF (lzu) zu = sum(work1(:,4)*cosmu(:))
259 IF (.not.lasym)
GOTO 1000
261 rmnsc = wmins*rzl_array(jslo,:,:,rsc)
262 1 + wplus*rzl_array(jshi,:,:,rsc)
263 zmncc = wmins*rzl_array(jslo,:,:,zcc)
264 1 + wplus*rzl_array(jshi,:,:,zcc)
267 rmncs = wmins*rzl_array(jslo,:,:,rcs)
268 1 + wplus*rzl_array(jshi,:,:,rcs)
269 zmnss = wmins*rzl_array(jslo,:,:,zss)
270 1 + wplus*rzl_array(jshi,:,:,zss)
280 work1(m,1) = sum(rmnsc(:,m)*cosnv(:))
281 work1(m,2) = sum(zmncc(:,m)*cosnv(:))
282 IF (lru) work1(m,3) = m*work1(m,1)
283 IF (lzu) work1(m,4) =-m*work1(m,2)
286 IF (lrv) work1(m,5) = sum(rmnsc(:,m)*sinnvn(:))
287 IF (lzv) work1(m,6) = sum(zmncc(:,m)*sinnvn(:))
288 work1(m,7) = sum(rmncs(:,m)*sinnv(:))
289 work1(m,8) = sum(zmnss(:,m)*sinnv(:))
290 IF (lru) work1(m,9) =-m*work1(m,7)
291 IF (lzu) work1(m,10) = m*work1(m,8)
292 IF (lrv) work1(m,11) = sum(rmncs(:,m)*cosnvn(:))
293 IF (lzv) work1(m,12) = sum(zmnss(:,m)*cosnvn(:))
302 r11 = r11 + sum(work1(:,1)*sinmu(:) + work1(:,7)*cosmu(:))
303 z11 = z11 + sum(work1(:,2)*cosmu(:) + work1(:,8)*sinmu(:))
305 1 sum(work1(:,3)*cosmu(:) + work1(:,9)*sinmu(:))
307 1 sum(work1(:,4)*sinmu(:) + work1(:,10)*cosmu(:))
309 1 sum(work1(:,5)*sinmu(:) + work1(:,11)*cosmu(:))
311 1 sum(work1(:,6)*cosmu(:) + work1(:,12)*sinmu(:))
313 r11 = r11 + sum(work1(:,1)*sinmu(:))
314 z11 = z11 + sum(work1(:,2)*cosmu(:))
315 IF (lru) ru = ru + sum(work1(:,3)*cosmu(:))
316 IF (lzu) zu = zu + sum(work1(:,4)*sinmu(:))
321 r_cyl(1) = r11; r_cyl(3) = z11
323 END SUBROUTINE flx2cyl
329 SUBROUTINE cyl2flx(rzl_in, r_cyl, c_flx, ns_in, ntor_in, mpol_in,
330 1 ntmax_in, lthreed_in, lasym_in, info, nfe, fmin,
331 1 mscale, nscale, ru, zu, rv, zv)
336 INTEGER,
INTENT(out) :: info, nfe
337 INTEGER,
INTENT(in) :: ns_in, ntor_in, mpol_in, ntmax_in
338 REAL(rprec),
INTENT(in) :: r_cyl(3)
339 REAL(rprec),
INTENT(inout) :: c_flx(3)
340 REAL(rprec),
INTENT(in),
TARGET ::
341 1 rzl_in(ns_in,0:ntor_in,0:mpol_in-1,3*ntmax_in)
342 REAL(rprec),
TARGET,
OPTIONAL ::
343 1 mscale(0:mpol_in-1), nscale(0:ntor_in)
344 REAL(rprec),
INTENT(out),
OPTIONAL :: ru, zu, rv, zv
345 REAL(rprec),
INTENT(out) :: fmin
346 LOGICAL,
INTENT(in) :: lthreed_in, lasym_in
350 INTEGER,
PARAMETER :: nvar = 2
351 real(rprec),
PARAMETER :: ftol = 1.e-16_dp
355 REAL(rprec) :: xc_opt(nvar), r_cyl_out(3), fmin0
356 INTEGER :: iflag, itry, nfe_out
361 INTEGER :: ns_loc, ntmax_loc, mpol_loc, ntor_loc
362 REAL(rprec) :: r_target, phi_target, z_target, fnorm
363 REAL(rprec),
POINTER :: rzl_array(:,:,:,:)
364 REAL(rprec),
POINTER :: mscale_loc(:), nscale_loc(:)
365 LOGICAL :: lthreed_loc, lasym_loc, lscale
404 lthreed_loc = lthreed_in; lasym_loc = lasym_in
405 mpol_loc = mpol_in; ntor_loc = ntor_in
406 ns_loc = ns_in; ntmax_loc = ntmax_in
407 lscale =
PRESENT(mscale)
409 mscale_loc => mscale; nscale_loc => nscale
411 r_target = r_cyl(1); phi_target = r_cyl(2); z_target = r_cyl(3)
414 xc_opt(1) = c_flx(1); xc_opt(2) = c_flx(2)
417 IF (c_flx(1) .eq. zero) xc_opt(1) = one/(ns_loc-1)
419 fnorm = r_target**2 + z_target**2
420 IF (fnorm .lt. epsilon(fnorm)) fnorm = 1
429 CALL newt2d(xc_opt, fmin, ftol, nfe_out, nvar, info)
432 IF (fmin.le.ftol .or. info.eq.-3)
EXIT
436 IF (fmin .gt. 1.e-3*fmin0)
THEN
437 xc_opt(2) = xc_opt(2) + twopi/20
439 xc_opt(2) = xc_opt(2) + twopi/40
442 fmin0 = min(fmin, fmin0)
447 c_flx(1) = xc_opt(1); c_flx(2) = xc_opt(2); c_flx(3) = phi_target
450 c_flx(2) = mod(c_flx(2), twopi)
451 DO WHILE (c_flx(2) .lt. zero)
452 c_flx(2) = c_flx(2) + twopi
465 IF (
PRESENT(ru) .or.
PRESENT(zu) .or.
PRESENT(rv) .or.
467 IF (info .eq. 0)
THEN
468 CALL flx2cyl(rzl_in, c_flx, r_cyl_out, ns_loc, ntor_loc,
469 & mpol_loc, ntmax_loc, lthreed_loc, lasym_loc,
470 & iflag, mscale=mscale, nscale=nscale,
471 & ru=ru, zu=zu, rv=rv, zv=zv)
473 IF (
PRESENT(ru)) ru = 0
474 IF (
PRESENT(zu)) zu = 0
475 IF (
PRESENT(rv)) rv = 0
476 IF (
PRESENT(zv)) zv = 0
484 SUBROUTINE newt2d(xc_opt, fmin, ftol, nfe, nvar, iflag)
489 INTEGER,
INTENT(in) :: nvar
490 INTEGER,
INTENT(out) :: nfe, iflag
491 REAL(rprec),
INTENT(inout) :: xc_opt(nvar)
492 REAL(rprec),
INTENT(in) :: ftol
493 REAL(rprec),
INTENT(out) :: fmin
497 INTEGER,
PARAMETER :: niter = 50
498 integer :: ieval, isgt1
499 REAL(rprec) :: c_flx(3), r_cyl_out(3), fvec(nvar), sflux,
500 1 uflux, eps0, eps, epu, xc_min(2), factor
501 REAL(rprec) :: x0(3), xs(3), xu(3), dels, delu, tau, fmin0,
502 1 ru1, zu1, edge_value, snew
530 eps0 = sqrt(epsilon(eps))
533 c_flx(3) = phi_target
537 edge_value = one + one/(ns_loc-1)
543 sflux = max(xc_opt(1), zero)
546 c_flx(1) = sflux; c_flx(2) = uflux
549 CALL get_flxcoord(x0, c_flx, ru=ru1, zu=zu1)
550 xu(1) = ru1; xu(3) = zu1
554 IF (sflux .ge. 1000*eps0)
THEN
562 IF (sflux .ge. 1-eps) eps = -eps
563 c_flx(1) = sflux + eps
564 CALL get_flxcoord(r_cyl_out, c_flx)
565 xs = (r_cyl_out - x0)/eps
568 x0(1) = x0(1) - r_target
569 x0(3) = x0(3) - z_target
570 fmin = (x0(1)**2 + x0(3)**2)*fnorm
572 IF (fmin .gt. fmin0)
THEN
573 factor = (2*factor)/3
576 IF (6*factor .lt. one)
THEN
577 dels =-(x0(1)*xs(1) + x0(3)*xs(3))/(xs(1)**2 + xs(3)**2)
578 delu =-(x0(1)*xu(1) + x0(3)*xu(3))/(xu(1)**2 + xu(3)**2)
586 tau = xu(1)*xs(3) - xu(3)*xs(1)
587 IF (abs(tau) .le. abs(eps)*r_target**2)
THEN
591 dels = ( x0(1)*xu(3) - x0(3)*xu(1))/tau
592 delu = (-x0(1)*xs(3) + x0(3)*xs(1))/tau
593 IF (fmin .gt. 1.e-3_dp)
THEN
594 dels = dels/2; delu = delu/2
599 IF (fmin .le. ftol)
EXIT
601 IF (abs(dels) .gt. one) dels = sign(one, dels)
604 snew = xc_opt(1) + dels*factor
605 IF (snew .lt. zero)
THEN
607 xc_opt(2) = xc_opt(2) + twopi/2
614 xc_opt(2) = xc_opt(2) + delu*factor
616 IF (xc_opt(1) .gt. edge_value)
THEN
618 IF (xc_opt(1) .gt. 2._dp) isgt1 = isgt1+1
619 IF (isgt1 .gt. 5)
EXIT
627 ELSE IF (xc_min(1) .gt. edge_value)
THEN
629 ELSE IF (fmin0 .le. ftol)
THEN
637 xc_opt(2) = mod(xc_opt(2), twopi)
639 END SUBROUTINE newt2d
643 SUBROUTINE get_flxcoord(x1, c_flx, ru, zu)
647 REAL(rprec),
INTENT(out) :: x1(3)
648 REAL(rprec),
INTENT(in) :: c_flx(3)
649 REAL(rprec),
INTENT(out),
OPTIONAL :: ru, zu
656 CALL flx2cyl(rzl_array, c_flx, x1, ns_loc, ntor_loc, mpol_loc,
657 1 ntmax_loc, lthreed_loc, lasym_loc, iflag,
658 2 mscale=mscale_loc, nscale=nscale_loc, ru=ru, zu=zu)
660 CALL flx2cyl(rzl_array, c_flx, x1, ns_loc, ntor_loc, mpol_loc,
661 1 ntmax_loc, lthreed_loc, lasym_loc, iflag,
665 END SUBROUTINE get_flxcoord
667 END SUBROUTINE cyl2flx