V3FIT
track_mod.f90
1 MODULE track_mod
2 !-------------------------------------------------------------------------------
3 !TRACK-TRACKs line segments through a torus
4 !
5 !TRACK_MOD is an F90 module of routines to obtain plasma information along a
6 ! segmented path
7 !
8 !References:
9 !
10 ! Attenberger, Houlberg, Hirshman, J Comp Phys 72 (1987) 435
11 ! W.A.Houlberg, P.I.Strand, D.McCune 4/2002
12 ! W.A.Houlberg, F90 free format 8/2004
13 !
14 !Contains PUBLIC routines:
15 !
16 ! TRACK -intersections of a segmented path with flux surfaces
17 ! TRACK_D -coordinates of point a distance d along segment
18 ! TRACK_G -Cartesian gradients along segment
19 !
20 !Comments:
21 !
22 ! Version 2.0 is a major upgrade of a previous package that has been in wide
23 ! use for about 15 years for RF, pellet, diagnostic and other tokamak and
24 ! stellarator applications.
25 !
26 ! The revisions include a faster and more robust algorithm for the
27 ! intersections.
28 !
29 ! This module is designed to map segmented straight-line trajectories in real
30 ! space to flux coordinates in order to find intersections with a
31 ! user-specified set of flux surfaces (TRACK).
32 !
33 ! Two auxiliary routines (TRACK_D and TRACK_G) are also accessible to the user
34 ! to find other information about a segment.
35 !
36 ! TRACK is meant to be used in transport analysis codes for mapping RF ray,
37 ! NBI, pellet, diagnostic or other trajectories to plasma coordinates for
38 ! plasma-neutral, plasma-wave, or other interactions that require information
39 ! from different coordinate systems.
40 !
41 ! The general tasks the TRACK modue performs are inherently 3D, so no
42 ! axisymmetric assumptions are made - that possibe restriction is left to the
43 ! MHD equilibrium interface module that actually performs the coordinate
44 ! conversions.
45 !
46 ! Two MHD equilibrium interface options are available:
47 ! AJAX -2D or 3D inverse coordinate representations implemented
48 ! by direct calls to AJAX routines
49 ! XPLASMA -2D inverse coordinate or Psi(R,Z) representations
50 ! implemented by calls to a module that maps AJAX calls to
51 ! XPLASMA calls (AJAX_XPLASMA_MOD)
52 !
53 ! The modernization of the code structure into an F90 module takes advantage of
54 ! some of the more attractive features of F90:
55 ! -use of KIND for precision declarations
56 ! -optional arguments for I/O
57 ! -generic names for all intrinsic functions
58 ! -compilation using either free or fixed form
59 ! -no common blocks or other deprecated Fortran features
60 ! -dynamic and automatic alocation of variables
61 ! -array syntax for vector operations
62 !-------------------------------------------------------------------------------
63 USE spec_kind_mod
64 
65 !Uncomment the appropriate MHD equilibrium interface, comment out the other
66 !For XPLASMA use mapping from XPLASMA to AJAX calls
67 !USE AJAX_XPLASMA_MOD
68 
69 !For AJAX use direct calls
70 USE ajax_mod
71 
72 IMPLICIT NONE
73 
74 REAL(kind=rspec), PARAMETER :: one=1, zero=0
75 
76 CONTAINS
77 
78 SUBROUTINE track(n_rho,rho,n_seg,r_seg, &
79  n_int,irho_int,s_int,iflag,message, &
80  K_SEG,IZONE_INT,RFLX_INT,RCYL_INT,RCAR_INT,SDOTB_INT, &
81  SDOTPHI_INT)
82 !-------------------------------------------------------------------------------
83 !TRACK finds the intersections of a segmented path with a set of flux surfaces
84 !
85 !References:
86 ! Attenberger, Houlberg, Hirshman, J Comp Phys 72 (1987) 435
87 ! W.A.Houlberg, P.I.Strand, D.McCune 4/2002
88 ! W.A.Houlberg, F90 free format 8/2004
89 !
90 !Comments:
91 ! It is assumed that rho(i) is increasing with i so that the i'th volume
92 ! element is bounded by rho(i) and rho(i+1)
93 ! The region outside rho(n_rho) is volume element n_rho
94 !-------------------------------------------------------------------------------
95 
96 !Declaration of input variables
97 INTEGER, INTENT(IN) :: &
98  n_rho, & !number of radial nodes [-]
99  n_seg !number of points defining the segments = segments+1 [-]
100 
101 REAL(KIND=rspec), INTENT(IN) :: &
102  rho(:), & !radial nodes [-]
103  r_seg(:,:) !coordinates defining the segments (see K_SEG)
104 
105 !Declaration of optional input variables
106 INTEGER, INTENT(IN), OPTIONAL :: &
107  K_SEG !option for coordinates specifying the segments [-]
108  !=1 flux coordinates (rho,theta,zeta)
109  !=2 Cartesian coordinates (x,y,z)
110  !=else default cylindrical coordinates (R,phi,Z)
111 
112 !Declaration of output variables
113 CHARACTER(LEN=1024), INTENT(OUT) :: &
114  message !warning or error message [character]
115 
116 INTEGER, INTENT(OUT) :: &
117  iflag, & !error and warning flag [-]
118  !=-1 warning
119  !=0 none
120  !=1 error
121  n_int, & !number of nodes + intersections [-]
122  irho_int(:) !radial nodes of intersections [-]
123  !=0 if point is a node of the ray
124 
125 REAL(KIND=rspec), INTENT(OUT) :: &
126  s_int(:) !length along path to intersections [m]
127 
128 !Declaration of optional output variables
129 INTEGER, INTENT(OUT), OPTIONAL :: &
130  IZONE_INT(:) !zone being entered [-]
131 
132 REAL(KIND=rspec), INTENT(OUT), OPTIONAL :: &
133  rflx_int(:,:), & !flux coordinates [rho,rad,rad]
134  rcyl_int(:,:), & !cylindrical coordinates [m,rad,m]
135  rcar_int(:,:), & !Cartesian coordinates [m,m,m]
136  sdotb_int(:), & !cos of angle between path and B [-]
137  sdotphi_int(:) !cos of angle between path and phi [-]
138 
139 !-------------------------------------------------------------------------------
140 !Declaration of local variables
141 INTEGER :: &
142  i,ilo,ihit,k,mhalf
143 
144 REAL(KIND=rspec) :: &
145  d0,d1,dseg,drds0,drds1,drds1f,drds_min,drho_min,ds,ds_max,ds_min,r000, &
146  rhomax,sseg,tau,tol,tol1
147 
148 REAL(KIND=rspec) :: &
149  r_car0(1:3),r_cyl0(1:3),r_flx0(1:3),g_cyl0(1:6), &
150  r_car1(1:3),r_cyl1(1:3),r_flx1(1:3),g_cyl1(1:6), &
151  r_cars(1:3),g_cars(1:3),r_cylc(1:3,1:n_seg),b_car(1:3)
152 
153 !-------------------------------------------------------------------------------
154 !Initialization
155 !-------------------------------------------------------------------------------
156 !General tolerence for min step sizes, resolution, convergence
157 tol=1.0e-4
158 tol1=1.0+tol
159 
160 !Get rho and characterstic major radius
161 iflag=0
162 message=''
163 CALL ajax_globals(iflag,message, &
164  rhomax=rhomax, &
165  r000=r000)
166 
167 !Check messages
168 IF(iflag /= 0) THEN
169 
170  message='TRACK(1) ' // message
171  IF(iflag > 0) GOTO 9999
172 
173 ENDIF
174 
175 !Set limits on step sizes relative to scale of torus
176 ds_min=tol*r000
177 ds_max=0.02*r000
178 drho_min=tol*rhomax
179 drds_min=drho_min/ds_max
180 
181 !Initialize local variables
182 drds0=0
183 drds1=0
184 
185 !Initialize output arrays
186 irho_int(:)=0
187 s_int(:)=0
188 IF(PRESENT(izone_int)) izone_int(:)=0
189 IF(PRESENT(rflx_int)) rflx_int(1:3,:)=0
190 IF(PRESENT(rcyl_int)) rcyl_int(1:3,:)=0
191 IF(PRESENT(rcar_int)) rcar_int(1:3,:)=0
192 IF(PRESENT(sdotb_int)) sdotb_int(:)=0
193 IF(PRESENT(sdotphi_int)) sdotphi_int(:)=0
194 
195 !Set option for coordinates defining segments
196 IF(PRESENT(k_seg)) THEN
197 
198  k=k_seg
199 
200 ELSE
201 
202  k=0
203 
204 ENDIF
205 
206 !Set cylindrical coordinates defining segments
207 IF(k == 1) THEN
208 
209  !Flux to cylindrical conversion
210  DO i=1,n_seg !Over chord points
211 
212  iflag=0
213  message=''
214  CALL ajax_flx2cyl(r_seg(1:3,i), &
215  r_cylc(1:3,i),iflag,message)
216 
217  !Check messages
218  IF(iflag /= 0) THEN
219 
220  message='TRACK(2) ' // message
221  IF(iflag > 0) GOTO 9999
222 
223  ENDIF
224 
225  ENDDO !Over chord points
226 
227 ELSEIF(k == 2) THEN
228 
229  !Cartesian to cylindrical conversion
230  DO i=1,n_seg !Over chord points
231 
232  CALL ajax_car2cyl(r_seg(1:3,i),r_cylc(1:3,i))
233 
234  ENDDO !Over chord points
235 
236 ELSE
237 
238  !Default cylindrical
239  r_cylc(1:3,1:n_seg)=r_seg(1:3,1:n_seg)
240 
241 ENDIF
242 
243 !-------------------------------------------------------------------------------
244 !Set starting point of first segment
245 !-------------------------------------------------------------------------------
246 !Distance along entire path
247 sseg=0
248 
249 !Cylindrical coordinates
250 r_cyl0(:)=r_cylc(:,1)
251 
252 !Cartesian coordinates
253 CALL ajax_cyl2car(r_cyl0,r_car0)
254 
255 !Flux coordinates
256 r_flx0(1)=rhomax
257 r_flx0(2)=atan2(-r_cyl0(3),r_cyl0(1)-r000)
258 r_flx0(3)=r_cyl0(2)
259 iflag=0
260 message=''
261 CALL ajax_cyl2flx(r_cyl0, &
262  r_flx0,iflag,message, &
263  g_cyl=g_cyl0)
264 
265 !Check messages
266 IF(iflag > 0) THEN
267 
268  !Assume starting point is outside plasma with bad Jacobian and proceed
269  r_flx0(1)=1.2*rhomax
270  r_flx0(2)=atan2(-r_cyl0(3),r_cyl0(1)-r000)
271  r_flx0(3)=r_cyl0(2)
272  g_cyl0(:)=(/one,zero,zero,zero,-1.2*rhomax,zero/)
273  iflag=0
274  message=''
275 
276 ENDIF
277 
278 !Find nearest flux surface between starting point and axis (may be on surface)
279 IF(r_flx0(1) > rho(n_rho)+drho_min) THEN
280 
281  !Outside largest surface of interest to user
282  ilo=n_rho
283  irho_int(1)=0
284 
285 ELSE
286 
287  !Find radial starting radial zone or intersection
288  loop_i: DO i=n_rho,1,-1 !Over radial nodes
289 
290  ilo=i
291 
292  IF(abs(r_flx0(1)-rho(ilo)) <= drho_min) THEN
293 
294  !On a surface
295  irho_int(1)=ilo
296  EXIT loop_i
297 
298  ELSEIF(r_flx0(1) >= rho(ilo)) THEN
299 
300  !ilo is the lower surface
301  irho_int(1)=0
302  EXIT loop_i
303 
304  ENDIF
305 
306  ENDDO loop_i !Over radial nodes
307 
308 ENDIF
309 
310 !Output starting point information
311 n_int=1
312 s_int(n_int)=0
313 IF(PRESENT(izone_int)) izone_int(n_int)=ilo
314 IF(PRESENT(rflx_int)) rflx_int(1:3,n_int)=r_flx0(1:3)
315 IF(PRESENT(rcyl_int)) rcyl_int(1:3,n_int)=r_cyl0(1:3)
316 IF(PRESENT(rcar_int)) rcar_int(1:3,n_int)=r_car0(1:3)
317 
318 CALL track_g(r_cylc(1,1),r_cylc(1,2),dseg,g_cars)
319 
320 IF(PRESENT(sdotb_int)) THEN
321 
322  iflag=0
323  message=''
324  b_car(:)=0
325  CALL ajax_b(r_flx0,r_cyl0,g_cyl0, &
326  iflag,message, &
327  b_car=b_car)
328 
329  !Check messages
330  IF(iflag > 0) THEN
331 
332  message='TRACK(3)/ERROR:failed to get B components'
333  GOTO 9999
334 
335  ENDIF
336 
337  IF(sum(b_car(1:3)) /= 0.0) THEN
338 
339  sdotb_int(n_int)=sum(g_cars(1:3)*b_car(1:3)) &
340  /sqrt(b_car(1)**2+b_car(2)**2+b_car(3)**2)
341 
342  ENDIF
343 
344 ENDIF
345 
346 IF(PRESENT(sdotphi_int)) THEN
347 
348  sdotphi_int(n_int)=-g_cars(1)*sin(r_cyl0(2))+g_cars(2)*cos(r_cyl0(2))
349 
350 ENDIF
351 
352 !Initialize flux coordinates at 1
353 r_flx1(:)=r_flx0(:)
354 
355 !-------------------------------------------------------------------------------
356 !Loop over segments
357 !-------------------------------------------------------------------------------
358 DO i=1,n_seg-1 !Over segments
359 
360  !Set Cartesian gradients of segment
361  CALL track_g(r_cylc(1,i),r_cylc(1,i+1),dseg,g_cars)
362 
363  !Cartesian coordinates of segment start
364  CALL ajax_cyl2car(r_cylc(:,i),r_cars)
365 
366  !Distance from start of segment
367  d0=0
368 
369  !Cylindrical coordinates and other information
370  iflag=0
371  message=''
372  CALL track_d(d0,r_cars,g_cars, &
373  r_flx0, &
374  r_car0,r_cyl0,drds0,g_cyl0,tau,iflag,message)
375 
376  !Make sure gradient exceeds minimum
377  IF(abs(drds0) <= drds_min) drds0=sign(one,drds0)*tol1*drds_min
378 
379  !Check messages
380  IF(iflag > 0 .OR. &
381  tau < 0.0) THEN
382 
383  !Assume d0 is outside plasma with bad Jacobian and proceed
384  r_flx0(:)=(/1.1*rhomax,zero,zero/)
385  drds0=0
386  iflag=0
387  message=''
388 
389  ENDIF
390 
391  !Set step halving counter, and initialize position of point 1
392  mhalf=0
393  d1=0
394 
395 !-------------------------------------------------------------------------------
396 !Step along segment: 0 is beginning of each step, 1 is end of step
397 !-------------------------------------------------------------------------------
398  loop_k: DO k=1,1000 !Over steps in segment
399 
400  !Check limiting number of steps
401  IF(k == 1000) THEN
402 
403  !Too many steps
404  iflag=1
405  message='TRACK(4)/ERROR:too many steps'
406  GOTO 9999
407 
408  ENDIF
409 
410  !Check conditions for termination
411  IF(mhalf > 5) THEN
412 
413  !Too many step halvings
414  iflag=1
415  message='TRACK(5)/ERROR:too many step halvings'
416  GOTO 9999
417 
418  ENDIF
419 
420  IF(d1 > dseg-ds_min) EXIT loop_k
421 
422  !Set step size
423  !Set ds from drho/ds at 0
424  IF(abs(drds0) < drds_min) THEN
425 
426  !Weak slope, limit step size to maximum
427  ds=ds_max
428 
429  ELSEIF(drds0 < 0.0) THEN
430 
431  !Heading inward
432  IF(ilo == 0) THEN
433 
434  !Special case where innermost surface is not axis
435  ds=ds_max
436 
437  ELSE
438 
439  !Try to step past ilo
440  ds=1.2*abs((r_flx0(1)-rho(ilo))/drds0)
441 
442  ENDIF
443 
444  ELSE
445 
446  IF(ilo < n_rho) THEN
447 
448  !Heading outward, try to step past ilo+1
449  ds=1.2*abs((rho(ilo+1)-r_flx0(1))/drds0)
450 
451  ELSE
452 
453  !Outside plasma, limit step to maximum
454  ds=ds_max
455 
456  ENDIF
457 
458  ENDIF
459 
460  !Check restrictions on step size
461  ds=min(ds,dseg-d0,ds_max)
462  ds=ds/real(2**mhalf,rspec)
463  ds=max(ds,ds_min)
464 
465  !Set position in segment at end of step
466  d1=d0+ds
467 
468  !Get parameters at 1
469  r_flx1(:)=r_flx0(:)
470  iflag=0
471  message=''
472  CALL track_d(d1,r_cars,g_cars, &
473  r_flx1, &
474  r_car1,r_cyl1,drds1,g_cyl1,tau,iflag,message)
475 
476  !Make sure gradient exceeds minimum
477  IF(abs(drds1) <= drds_min) drds1=sign(one,drds1)*tol1*drds_min
478 
479  !Check messages
480  IF(iflag > 0 .OR. &
481  tau < 0.0) THEN
482 
483  !Assume d1 is outside plasma with bad Jacobian and proceed
484  r_flx1(:)=(/1.1*rhomax,zero,r_cyl1(2)/)
485  drds1=0
486  iflag=0
487  message=''
488 
489  ENDIF
490 
491 !-------------------------------------------------------------------------------
492 !Checks for intersections and tangencies depend on drho/ds at endpoints
493 !-------------------------------------------------------------------------------
494  IF(abs(drds0) < drds_min) THEN
495 
496 !-------------------------------------------------------------------------------
497 !Case I: Invalid soln at 0
498 
499  IF(abs(drds1) < drds_min) THEN
500 
501 !---------
502 !Case I.A: Invalid solns at 0 and 1
503  !Proceed with next step
504  ihit=0
505  mhalf=0
506 
507  ELSE
508 
509  IF(r_flx1(1) > rho(n_rho)) THEN
510 
511 !---------
512 !Case I.B&C.1: Outside target surface
513  !Proceed with next step
514  ihit=0
515  mhalf=0
516 
517  ELSE
518 
519 !---------
520 !Case I.B&C.2: Crossed target surface
521  !Reduce step size to get good data point outside plasma
522  ihit=0
523  mhalf=1
524 
525  ENDIF
526 
527  ENDIF
528 
529  ELSEIF(drds0 > 0.0) THEN
530 
531 !-------------------------------------------------------------------------------
532 !Case II: Valid outward soln at 0
533 
534  IF(abs(drds1) < drds_min) THEN
535 
536 !---------
537 !Case II.A: Valid outward soln at 0, invalid soln at 1
538 
539  IF(ilo == n_rho) THEN
540 
541 !Case II.A.1: Outside plasma
542  !Proceed with next step
543  ihit=0
544  mhalf=0
545 
546  ELSE
547 
548 !Case II.A.2: Inside plasma
549  !Reduce step to find intersection
550  ihit=0
551  mhalf=mhalf+1
552 
553  ENDIF
554 
555  ELSEIF(drds1 > 0.0) THEN
556 
557 !---------
558 !Case II.B: Valid outward soln at 0, valid outward soln at 1
559  IF(ilo == n_rho) THEN
560 
561 !Case II.B.1: Outside plasma
562  !Proceed with next step
563  ihit=0
564  mhalf=0
565 
566  ELSEIF(abs(r_flx1(1)-rho(ilo+1)) < drho_min) THEN
567 
568 !Case II.B.2: Hit surface ilo+1
569  !Record hit and proceed with next step
570  ihit=ilo+1
571  mhalf=0
572  ilo=ilo+1
573  d1=d1+max(ds_min,(rho(ihit)-r_flx1(1))/drds1)
574  ds=d1-d0
575 
576  ELSEIF(ilo < n_rho-1 .AND. &
577  r_flx1(1) > rho(ilo+2)) THEN
578 
579 !Case II.B.3: Crossed more than one surface
580  !Halve step
581  ihit=0
582  mhalf=mhalf+1
583 
584  ELSEIF(r_flx1(1) > rho(ilo+1)) THEN
585 
586 !Case II.B.4: Crossed surface ilo+1
587  !Get intersection, record hit, and proceed with next step
588  ihit=ilo+1
589  mhalf=0
590  ilo=ilo+1
591  d1=d0+max(ds_min,(ds*(rho(ihit)-r_flx0(1)) &
592  /(r_flx1(1)-r_flx0(1))))
593  ds=d1-d0
594 
595  ELSE
596 !Case II.B.5: Crossed no surfaces
597  !Proceed with next step
598  ihit=0
599  mhalf=0
600 
601  ENDIF
602 
603  ELSE
604 
605 !---------
606 !Case II.C: Valid outward soln at 0, valid inward soln at 1, passed tangency
607  !Get position of tangency
608  d1=max(d0+ds_min,(d1*drds0-d0*drds1)/(drds0-drds1))
609  ds=d1-d0
610 
611  !Get parameters at tangency
612  r_flx1(:)=r_flx0(:)
613  iflag=0
614  message=''
615  CALL track_d(d1,r_cars,g_cars, &
616  r_flx1, &
617  r_car1,r_cyl1,drds1,g_cyl1,tau,iflag,message)
618 
619  !Check messages
620  IF(iflag > 0 .OR. &
621  tau < 0.0) THEN
622 
623  message='TRACK(6) ' // message
624  GOTO 9999
625 
626  ENDIF
627 
628  !Force gradient to minimum beyond tangency
629  drds1=-drds_min
630 
631  !Check for intersection at tangency
632  IF(ilo == n_rho) THEN
633 
634 !Case II.C.1: Tangency outside outermost surface
635  !Proceed with next step
636  ihit=0
637  mhalf=0
638 
639  ELSEIF(abs(r_flx1(1)-rho(ilo+1)) < drho_min) THEN
640 
641 !Case II.C.2: Tangent to surface ilo+1 (degenerate case)
642  !Record hit and proceed with next step
643  ihit=ilo+1
644  mhalf=0
645 
646  ELSEIF(r_flx1(1) > rho(ilo+1)) THEN
647 
648 !Case II.C.3: Crossed surface ilo+1 before tangency
649  !Get intersection and proceed with next step
650  ihit=ilo+1
651  mhalf=0
652  ilo=ilo+1
653  d1=d0+max(ds_min,(ds*(rho(ihit)-r_flx0(1)) &
654  /(r_flx1(1)-r_flx0(1))))
655  ds=d1-d0
656 
657  !Unforce gradient and proceed outward
658  drds1=2*drds_min
659 
660  ELSE
661 
662 !Case II.C.4: Crossed no surface
663  !Proceed with next step
664  ihit=0
665  mhalf=0
666 
667  ENDIF
668 
669  ENDIF
670 
671  ELSE
672 
673 !-------------------------------------------------------------------------------
674 !Case III: Valid inward soln at 0
675 
676  IF(abs(drds1) < drds_min) THEN
677 
678 !---------
679 !Case III.A: Valid inward soln at 0, invalid soln at 1
680  IF(ilo == n_rho) THEN
681 
682 !Case III.A.1: Outside plasma
683  !Proceed with next step
684  ihit=0
685  mhalf=0
686 
687  ELSE
688 
689 !Case III.A.2: Inside plasma
690  !Failure
691  iflag=1
692  message='TRACK(7)/ERROR:conversion failure in plasma'
693  GOTO 9999
694 
695  ENDIF
696 
697  ELSEIF(drds1 > 0) THEN
698 
699 !---------
700 !Case III.B: Valid inward soln at 0, valid outward soln at 1, passed tangency
701  !Get position of tangency
702  d1=max(d0+ds_min,(d1*drds0-d0*drds1)/(drds0-drds1))
703  ds=d1-d0
704 
705  !Get parameters at tangency
706  r_flx1(:)=r_flx0(:)
707  iflag=0
708  message=''
709  CALL track_d(d1,r_cars,g_cars, &
710  r_flx1, &
711  r_car1,r_cyl1,drds1,g_cyl1,tau,iflag,message)
712 
713  !Force gradient to minimum beyond tangency
714  drds1=drds_min
715 
716  !Check messages
717  IF(iflag > 0 .OR. &
718  tau < 0.0) THEN
719 
720  message='TRACK(8) ' // message
721  GOTO 9999
722 
723  ENDIF
724 
725  !Check for intersection at tangency
726  IF(ilo == 0) THEN
727 
728 !Case III.B.1: Tangency inside innermost surface
729  !Proceed with next step
730  ihit=0
731  mhalf=0
732 
733  ELSEIF(rho(ilo) < drho_min) THEN
734 
735 !Case III.B.2: Don't look for tangency near axis
736  !Proceed with next step
737  ihit=0
738  mhalf=0
739 
740  ELSEIF(abs(r_flx1(1)-rho(ilo)) < drho_min) THEN
741 
742 !Case III.B.3: Tangent to surface ilo (degenerate case)
743  !Record hit and proceed with next step
744  ihit=ilo
745  mhalf=0
746 
747  ELSEIF(r_flx1(1) < rho(ilo)) THEN
748 
749 !Case III.B.4: Crossed surface ilo before tangency
750  !Get and record intersection and proceed with next step
751  ihit=ilo
752  mhalf=0
753  ilo=ilo-1
754  d1=d0+max(ds_min,(ds*(rho(ihit)-r_flx0(1)) &
755  /(r_flx1(1)-r_flx0(1))))
756  ds=d1-d0
757 
758  !Unforce gradient and proceed inward
759  drds1=-2*drds_min
760 
761  ELSE
762 
763 !Case III.B.5: Crossed no surface
764  !Proceed with next step
765  ihit=0
766  mhalf=0
767 
768  ENDIF
769 
770  ELSE
771 
772 !---------
773 !Case III.C: Valid inward soln at 0, valid inward soln at 1
774  IF(ilo == 0) THEN
775 
776 !Case III.C.1: Hasn't passed tangency to re-hit lowest surface
777  !Proceed with next step
778  ihit=0
779  mhalf=0
780 
781  ELSEIF(rho(ilo) < drho_min) THEN
782 
783 !Case III.C.2: Don't look for intersection near axis
784  !Proceed with next step
785  ihit=0
786  mhalf=0
787 
788  ELSEIF(abs(r_flx1(1)-rho(ilo)) < drho_min) THEN
789 
790 !Case III.C.3: Hit surface ilo
791  !Record hit and proceed with next step
792  ihit=ilo
793  mhalf=0
794  ilo=ilo-1
795  d1=d1+max(ds_min,(rho(ihit)-r_flx1(1))/drds1)
796  ds=d1-d0
797 
798  ELSEIF(ilo > 1 .AND. &
799  r_flx1(1) < rho(ilo-1)) THEN
800 
801 !Case III.C.4: Crossed more than one surface
802  !Halve step
803  ihit=0
804  mhalf=mhalf+1
805 
806  ELSEIF(r_flx1(1) < rho(ilo)) THEN
807 
808 !Case III.C.5: Crossed surface ilo
809  !Get intersection, record hit, and proceed with next step
810  ihit=ilo
811  mhalf=0
812  ilo=ilo-1
813  d1=d0+max(ds_min,(ds*(rho(ihit)-r_flx0(1)) &
814  /(r_flx1(1)-r_flx0(1))))
815  ds=d1-d0
816 
817  ELSE
818 
819 !Case III.C.6: Crossed no surface
820  !Proceed with next step
821  ihit=0
822  mhalf=0
823 
824  ENDIF
825 
826  ENDIF
827 
828  ENDIF
829 
830 !-------------------------------------------------------------------------------
831 !Step completed
832 !-------------------------------------------------------------------------------
833 
834  IF(ihit /= 0) THEN
835 
836  !Save in case of forced end gradients
837  drds1f=drds1
838 
839  !Record hit
840  iflag=0
841  message=''
842  CALL track_d(d1,r_cars,g_cars, &
843  r_flx1, &
844  r_car1,r_cyl1,drds1,g_cyl1,tau,iflag,message)
845 
846  IF(abs(drds1) <= tol1*drds_min) drds1=drds1f
847 
848  !Check messages
849  IF(iflag > 0 .OR. &
850  tau < 0.0) THEN
851 
852  message='TRACK(9) ' // message
853  GOTO 9999
854 
855  ENDIF
856 
857  n_int=n_int+1
858  r_flx1(1)=rho(ihit)
859  irho_int(n_int)=ihit
860  s_int(n_int)=sseg+d1
861  IF(PRESENT(izone_int) .AND. &
862  drds0 > 0.0) izone_int(n_int)=ihit
863  IF(PRESENT(izone_int) .AND. &
864  drds0 < 0.0) izone_int(n_int)=ihit-1
865  IF(PRESENT(rflx_int)) rflx_int(1:3,n_int)=r_flx1(1:3)
866  IF(PRESENT(rcyl_int)) rcyl_int(1:3,n_int)=r_cyl1(1:3)
867  IF(PRESENT(rcar_int)) rcar_int(1:3,n_int)=r_car1(1:3)
868 
869  IF(PRESENT(sdotb_int)) THEN
870 
871  iflag=0
872  message=''
873  b_car(:)=0
874  CALL ajax_b(r_flx1,r_cyl1,g_cyl1, &
875  iflag,message, &
876  b_car=b_car)
877 
878  !Check messages
879  IF(iflag > 0) THEN
880 
881  message='TRACK(10)/ERROR:failed to get B components'
882  GOTO 9999
883 
884  ENDIF
885 
886  sdotb_int(n_int)=sum(g_cars(1:3)*b_car(1:3)) &
887  /sqrt(b_car(1)**2+b_car(2)**2+b_car(3)**2)
888 
889  ENDIF
890 
891  IF(PRESENT(sdotphi_int)) THEN
892 
893  sdotphi_int(n_int)=-g_cars(1)*sin(r_cyl1(2))+g_cars(2) &
894  *cos(r_cyl1(2))
895 
896  ENDIF
897 
898  ihit=0
899 
900  ENDIF
901 
902  IF(mhalf == 0) THEN
903 
904  !Advance point 0 to point 1
905  r_flx0(1:3)=r_flx1(1:3)
906  r_cyl0(1:3)=r_cyl1(1:3)
907  r_car0(1:3)=r_car1(1:3)
908  g_cyl0(1:6)=g_cyl1(1:6)
909  drds0=drds1
910  d0=d1
911 
912  ENDIF
913 
914  ENDDO loop_k !Over steps in segment
915 
916 !-------------------------------------------------------------------------------
917 !Segment completed
918 !-------------------------------------------------------------------------------
919  sseg=sseg+dseg
920  n_int=n_int+1
921  irho_int(n_int)=0
922  s_int(n_int)=sseg
923  r_cyl1(1:3)=r_cylc(1:3,i+1)
924 
925  iflag=0
926  message=''
927  CALL ajax_cyl2flx(r_cyl1, &
928  r_flx1,iflag,message, &
929  g_cyl=g_cyl1)
930 
931  !Check messages
932  IF(iflag > 0) THEN
933 
934  !Assume d1 is outside plasma with bad Jacobian and proceed
935  r_flx1(:)=(/1.1*rhomax,zero,zero/)
936  drds1=0
937  iflag=-1
938  message='TRACK(11)/WARNING:end point assumed outside plasma'
939 
940  ENDIF
941 
942  CALL ajax_cyl2car(r_cyl1,r_car1)
943 
944  IF(PRESENT(izone_int)) izone_int(n_int)=ilo
945  IF(PRESENT(rflx_int)) rflx_int(1:3,n_int)=r_flx1(1:3)
946  IF(PRESENT(rcyl_int)) rcyl_int(1:3,n_int)=r_cyl1(1:3)
947  IF(PRESENT(rcar_int)) rcar_int(1:3,n_int)=r_car1(1:3)
948 
949  IF(PRESENT(sdotb_int)) THEN
950 
951  iflag=0
952  message=''
953  CALL ajax_b(r_flx1,r_cyl1,g_cyl1, &
954  iflag,message, &
955  b_car=b_car)
956 
957  !Check messages
958  IF(iflag > 0) THEN
959 
960  message='TRACK(12)/ERROR:failed to get B components'
961  GOTO 9999
962 
963  ENDIF
964 
965  sdotb_int(n_int)=sum(g_cars(1:3)*b_car(1:3)) &
966  /sqrt(b_car(1)**2+b_car(2)**2+b_car(3)**2)
967 
968  ENDIF
969 
970  IF(PRESENT(sdotphi_int)) THEN
971 
972  sdotphi_int(n_int)=-g_cars(1)*sin(r_cylc(2,i+1)) &
973  +g_cars(2)*cos(r_cylc(2,i+1))
974 
975  ENDIF
976 
977 ENDDO !Over segments
978 
979 !Make sure intersection distances are monotonic
980 DO k=1,n_int-1
981 
982  IF(s_int(k+1) <= s_int(k)) THEN
983 
984  iflag=1
985  message='TRACK(12)/ERROR:intercept lengths out of order.'
986  GOTO 9999
987 
988  ENDIF
989 
990 ENDDO
991 
992 !-------------------------------------------------------------------------------
993 !Cleanup and exit
994 !-------------------------------------------------------------------------------
995 9999 CONTINUE
996 
997 END SUBROUTINE track
998 
999 SUBROUTINE track_d(d,r_cars,g_cars, &
1000  r_flx, &
1001  r_car,r_cyl,drhods,g_cyl,tau,iflag,message)
1002 !-------------------------------------------------------------------------------
1003 !TRACK_D finds the Cartesian, cylindrical and flux coordinates of a point
1004 ! located a distance d along a line segment specified in Cartesian coordinates,
1005 ! and returns other relevant information for that point
1006 !
1007 !References:
1008 ! Attenberger, Houlberg, Hirshman J Comp Phys 72 (1987) 435
1009 ! W.A.Houlberg, P.I.Strand 11/2001
1010 ! W.A.Houlberg, F90 free format 8/2004
1011 !
1012 !Comments:
1013 ! For efficiency, r_flx should be a good initial guess on entry
1014 !-------------------------------------------------------------------------------
1015 
1016 !Declaration of input variables
1017 REAL(KIND=rspec), INTENT(IN) :: &
1018  d, & !distance along segment [m]
1019  r_cars(:), & !Cartesian coordinates of start of segment [m]
1020  g_cars(:) !unit length vector along segment [-]
1021 
1022 !Declaration of input/output variables
1023 REAL(KIND=rspec), INTENT(INOUT) :: &
1024  r_flx(:) !flux coordinates (rho,theta,zeta) [rho,rad,rad]
1025 
1026 !Declaration of output variables
1027 CHARACTER(LEN=1024), INTENT(OUT) :: &
1028  message !warning or error message [character]
1029 
1030 INTEGER, INTENT(OUT) :: &
1031  iflag !error and warning flag [-]
1032  !=-1 warning
1033  !=0 none
1034  !=1 error
1035 
1036 
1037 REAL(KIND=rspec), INTENT(OUT) :: &
1038  r_car(:), & !Cartesian coordinates (x,y,z) [m,m,m]
1039  r_cyl(:), & !cylindrical coordinates (R,phi,Z) [m,rad,m]
1040  g_cyl(:), & !R,Z derivatives
1041  !=(R_rho,R_theta,R_zeta,Z_rho,Z_theta,Z_zeta)
1042  ! [m/rho,m, m, m/rho,m, m ]
1043  drhods, & !drho/ds at position d along the chord [rho/m]
1044  tau !2-D Jacobian in phi=zeta=constant plane [m**2/rho]
1045 
1046 !-------------------------------------------------------------------------------
1047 !Declaration of local variables
1048 REAL(KIND=rspec) :: &
1049  rhop,rhor,rhoz
1050 
1051 !-------------------------------------------------------------------------------
1052 !Get coordinates
1053 !-------------------------------------------------------------------------------
1054 !Cartesian coordinates of end point
1055 r_car(1:3)=r_cars(1:3)+g_cars(1:3)*d
1056 
1057 !Cylindrical coordinates
1058 CALL ajax_car2cyl(r_car, &
1059  r_cyl)
1060 
1061 !Flux coordinates
1062 iflag=0
1063 message=''
1064 CALL ajax_cyl2flx(r_cyl, &
1065  r_flx,iflag,message, &
1066  g_cyl=g_cyl, &
1067  tau=tau)
1068 
1069 !Check messages
1070 IF(iflag >0) THEN
1071 
1072  message='TRACK_D ' // message
1073 
1074 ENDIF
1075 
1076 !drho/ds
1077 rhor=-g_cyl(5)/tau
1078 rhoz=g_cyl(2)/tau
1079 rhop=(g_cyl(3)*g_cyl(5)-g_cyl(2)*g_cyl(6))/tau
1080 drhods=(rhor*r_car(1)-rhop*r_car(2)/r_cyl(1))*g_cars(1)/r_cyl(1) &
1081  +(rhor*r_car(2)+rhop*r_car(1)/r_cyl(1))*g_cars(2)/r_cyl(1)+rhoz*g_cars(3)
1082 
1083 END SUBROUTINE track_d
1084 
1085 SUBROUTINE track_g(r_cyl0,r_cyl1, &
1086  d,g_car)
1087 !-------------------------------------------------------------------------------
1088 !TRACK_G evaluates the Cartesian gradients along a segment defined by two sets
1089 ! of cylindrical coordinates
1090 !
1091 !References:
1092 ! Attenberger, Houlberg, Hirshman, J Comp Phys 72 (1987) 435
1093 ! W.A.Houlberg 8/2001
1094 ! W.A.Houlberg, F90 free format 8/2004
1095 !-------------------------------------------------------------------------------
1096 
1097 !Declaration of input variables
1098 REAL(KIND=rspec), INTENT(IN) :: &
1099  r_cyl0(3), & !cylindrical coordinates at start of segment [m,rad,m]
1100  r_cyl1(3) !cylindrical coordinates at end of segment [m,rad,m]
1101 
1102 !Declaration of output variables
1103 REAL(KIND=rspec), INTENT(OUT) :: &
1104  d, & !length of segment [m]
1105  g_car(3) !unit length vector along segment [-]
1106 
1107 !-------------------------------------------------------------------------------
1108 !Declaration of local variables
1109 REAL(KIND=rspec) :: &
1110  r_car0(3),r_car1(3)
1111 
1112 !-------------------------------------------------------------------------------
1113 !Get gradients
1114 !-------------------------------------------------------------------------------
1115 !Cartesian coordinates for endpoints
1116 CALL ajax_cyl2car(r_cyl0,r_car0)
1117 CALL ajax_cyl2car(r_cyl1,r_car1)
1118 
1119 !Increments in Cartesian coordinates
1120 g_car(:)=r_car1(:)-r_car0(:)
1121 
1122 !Segment length
1123 d=sqrt((r_car1(1)-r_car0(1))**2+(r_car1(2)-r_car0(2))**2 &
1124  +(r_car1(3)-r_car0(3))**2)
1125 
1126 !Normalized gradients
1127 g_car(:)=g_car(:)/d
1128 
1129 END SUBROUTINE track_g
1130 
1131 END MODULE track_mod