74 REAL(kind=rspec),
PARAMETER :: one=1, zero=0
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, &
97 INTEGER,
INTENT(IN) :: &
101 REAL(KIND=rspec),
INTENT(IN) :: &
106 INTEGER,
INTENT(IN),
OPTIONAL :: &
113 CHARACTER(LEN=1024),
INTENT(OUT) :: &
116 INTEGER,
INTENT(OUT) :: &
125 REAL(KIND=rspec),
INTENT(OUT) :: &
129 INTEGER,
INTENT(OUT),
OPTIONAL :: &
132 REAL(KIND=rspec),
INTENT(OUT),
OPTIONAL :: &
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
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)
163 CALL ajax_globals(iflag,message, &
170 message=
'TRACK(1) ' // message
171 IF(iflag > 0)
GOTO 9999
179 drds_min=drho_min/ds_max
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
196 IF(
PRESENT(k_seg))
THEN
214 CALL ajax_flx2cyl(r_seg(1:3,i), &
215 r_cylc(1:3,i),iflag,message)
220 message=
'TRACK(2) ' // message
221 IF(iflag > 0)
GOTO 9999
232 CALL ajax_car2cyl(r_seg(1:3,i),r_cylc(1:3,i))
239 r_cylc(1:3,1:n_seg)=r_seg(1:3,1:n_seg)
250 r_cyl0(:)=r_cylc(:,1)
253 CALL ajax_cyl2car(r_cyl0,r_car0)
257 r_flx0(2)=atan2(-r_cyl0(3),r_cyl0(1)-r000)
261 CALL ajax_cyl2flx(r_cyl0, &
262 r_flx0,iflag,message, &
270 r_flx0(2)=atan2(-r_cyl0(3),r_cyl0(1)-r000)
272 g_cyl0(:)=(/one,zero,zero,zero,-1.2*rhomax,zero/)
279 IF(r_flx0(1) > rho(n_rho)+drho_min)
THEN
288 loop_i:
DO i=n_rho,1,-1
292 IF(abs(r_flx0(1)-rho(ilo)) <= drho_min)
THEN
298 ELSEIF(r_flx0(1) >= rho(ilo))
THEN
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)
318 CALL track_g(r_cylc(1,1),r_cylc(1,2),dseg,g_cars)
320 IF(
PRESENT(sdotb_int))
THEN
325 CALL ajax_b(r_flx0,r_cyl0,g_cyl0, &
332 message=
'TRACK(3)/ERROR:failed to get B components'
337 IF(sum(b_car(1:3)) /= 0.0)
THEN
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)
346 IF(
PRESENT(sdotphi_int))
THEN
348 sdotphi_int(n_int)=-g_cars(1)*sin(r_cyl0(2))+g_cars(2)*cos(r_cyl0(2))
361 CALL track_g(r_cylc(1,i),r_cylc(1,i+1),dseg,g_cars)
364 CALL ajax_cyl2car(r_cylc(:,i),r_cars)
372 CALL track_d(d0,r_cars,g_cars, &
374 r_car0,r_cyl0,drds0,g_cyl0,tau,iflag,message)
377 IF(abs(drds0) <= drds_min) drds0=sign(one,drds0)*tol1*drds_min
384 r_flx0(:)=(/1.1*rhomax,zero,zero/)
405 message=
'TRACK(4)/ERROR:too many steps'
415 message=
'TRACK(5)/ERROR:too many step halvings'
420 IF(d1 > dseg-ds_min)
EXIT loop_k
424 IF(abs(drds0) < drds_min)
THEN
429 ELSEIF(drds0 < 0.0)
THEN
440 ds=1.2*abs((r_flx0(1)-rho(ilo))/drds0)
449 ds=1.2*abs((rho(ilo+1)-r_flx0(1))/drds0)
461 ds=min(ds,dseg-d0,ds_max)
462 ds=ds/real(2**mhalf,rspec)
472 CALL track_d(d1,r_cars,g_cars, &
474 r_car1,r_cyl1,drds1,g_cyl1,tau,iflag,message)
477 IF(abs(drds1) <= drds_min) drds1=sign(one,drds1)*tol1*drds_min
484 r_flx1(:)=(/1.1*rhomax,zero,r_cyl1(2)/)
494 IF(abs(drds0) < drds_min)
THEN
499 IF(abs(drds1) < drds_min)
THEN
509 IF(r_flx1(1) > rho(n_rho))
THEN
529 ELSEIF(drds0 > 0.0)
THEN
534 IF(abs(drds1) < drds_min)
THEN
539 IF(ilo == n_rho)
THEN
555 ELSEIF(drds1 > 0.0)
THEN
559 IF(ilo == n_rho)
THEN
566 ELSEIF(abs(r_flx1(1)-rho(ilo+1)) < drho_min)
THEN
573 d1=d1+max(ds_min,(rho(ihit)-r_flx1(1))/drds1)
576 ELSEIF(ilo < n_rho-1 .AND. &
577 r_flx1(1) > rho(ilo+2))
THEN
584 ELSEIF(r_flx1(1) > rho(ilo+1))
THEN
591 d1=d0+max(ds_min,(ds*(rho(ihit)-r_flx0(1)) &
592 /(r_flx1(1)-r_flx0(1))))
608 d1=max(d0+ds_min,(d1*drds0-d0*drds1)/(drds0-drds1))
615 CALL track_d(d1,r_cars,g_cars, &
617 r_car1,r_cyl1,drds1,g_cyl1,tau,iflag,message)
623 message=
'TRACK(6) ' // message
632 IF(ilo == n_rho)
THEN
639 ELSEIF(abs(r_flx1(1)-rho(ilo+1)) < drho_min)
THEN
646 ELSEIF(r_flx1(1) > rho(ilo+1))
THEN
653 d1=d0+max(ds_min,(ds*(rho(ihit)-r_flx0(1)) &
654 /(r_flx1(1)-r_flx0(1))))
676 IF(abs(drds1) < drds_min)
THEN
680 IF(ilo == n_rho)
THEN
692 message=
'TRACK(7)/ERROR:conversion failure in plasma'
697 ELSEIF(drds1 > 0)
THEN
702 d1=max(d0+ds_min,(d1*drds0-d0*drds1)/(drds0-drds1))
709 CALL track_d(d1,r_cars,g_cars, &
711 r_car1,r_cyl1,drds1,g_cyl1,tau,iflag,message)
720 message=
'TRACK(8) ' // message
733 ELSEIF(rho(ilo) < drho_min)
THEN
740 ELSEIF(abs(r_flx1(1)-rho(ilo)) < drho_min)
THEN
747 ELSEIF(r_flx1(1) < rho(ilo))
THEN
754 d1=d0+max(ds_min,(ds*(rho(ihit)-r_flx0(1)) &
755 /(r_flx1(1)-r_flx0(1))))
781 ELSEIF(rho(ilo) < drho_min)
THEN
788 ELSEIF(abs(r_flx1(1)-rho(ilo)) < drho_min)
THEN
795 d1=d1+max(ds_min,(rho(ihit)-r_flx1(1))/drds1)
798 ELSEIF(ilo > 1 .AND. &
799 r_flx1(1) < rho(ilo-1))
THEN
806 ELSEIF(r_flx1(1) < rho(ilo))
THEN
813 d1=d0+max(ds_min,(ds*(rho(ihit)-r_flx0(1)) &
814 /(r_flx1(1)-r_flx0(1))))
842 CALL track_d(d1,r_cars,g_cars, &
844 r_car1,r_cyl1,drds1,g_cyl1,tau,iflag,message)
846 IF(abs(drds1) <= tol1*drds_min) drds1=drds1f
852 message=
'TRACK(9) ' // message
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)
869 IF(
PRESENT(sdotb_int))
THEN
874 CALL ajax_b(r_flx1,r_cyl1,g_cyl1, &
881 message=
'TRACK(10)/ERROR:failed to get B components'
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)
891 IF(
PRESENT(sdotphi_int))
THEN
893 sdotphi_int(n_int)=-g_cars(1)*sin(r_cyl1(2))+g_cars(2) &
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)
923 r_cyl1(1:3)=r_cylc(1:3,i+1)
927 CALL ajax_cyl2flx(r_cyl1, &
928 r_flx1,iflag,message, &
935 r_flx1(:)=(/1.1*rhomax,zero,zero/)
938 message=
'TRACK(11)/WARNING:end point assumed outside plasma'
942 CALL ajax_cyl2car(r_cyl1,r_car1)
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)
949 IF(
PRESENT(sdotb_int))
THEN
953 CALL ajax_b(r_flx1,r_cyl1,g_cyl1, &
960 message=
'TRACK(12)/ERROR:failed to get B components'
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)
970 IF(
PRESENT(sdotphi_int))
THEN
972 sdotphi_int(n_int)=-g_cars(1)*sin(r_cylc(2,i+1)) &
973 +g_cars(2)*cos(r_cylc(2,i+1))
982 IF(s_int(k+1) <= s_int(k))
THEN
985 message=
'TRACK(12)/ERROR:intercept lengths out of order.'
999 SUBROUTINE track_d(d,r_cars,g_cars, &
1001 r_car,r_cyl,drhods,g_cyl,tau,iflag,message)
1017 REAL(KIND=rspec),
INTENT(IN) :: &
1023 REAL(KIND=rspec),
INTENT(INOUT) :: &
1027 CHARACTER(LEN=1024),
INTENT(OUT) :: &
1030 INTEGER,
INTENT(OUT) :: &
1037 REAL(KIND=rspec),
INTENT(OUT) :: &
1048 REAL(KIND=rspec) :: &
1055 r_car(1:3)=r_cars(1:3)+g_cars(1:3)*d
1058 CALL ajax_car2cyl(r_car, &
1064 CALL ajax_cyl2flx(r_cyl, &
1065 r_flx,iflag,message, &
1072 message=
'TRACK_D ' // message
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)
1083 END SUBROUTINE track_d
1085 SUBROUTINE track_g(r_cyl0,r_cyl1, &
1098 REAL(KIND=rspec),
INTENT(IN) :: &
1103 REAL(KIND=rspec),
INTENT(OUT) :: &
1109 REAL(KIND=rspec) :: &
1116 CALL ajax_cyl2car(r_cyl0,r_car0)
1117 CALL ajax_cyl2car(r_cyl1,r_car1)
1120 g_car(:)=r_car1(:)-r_car0(:)
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)
1129 END SUBROUTINE track_g
1131 END MODULE track_mod