V3FIT
rt_hand_side_dphi.f
1  SUBROUTINE rt_hand_side_dphi(neq, phi, yin, fout) ! Defines field-line equations for LSODE integrator for type="toroidal"
2  USE stel_kinds ! The right hand sides are B^s/B^v and B^u/B^v.
3  USE stel_constants ! To get them, one simply uses:
4  USE bspline ! B^R =B*grad R = B^sRs + B^uRu + B^vRv
5  USE spline3d_fit_coefs ! B^z =B*grad Z = B^sZs + B^uZu + B^vZu
6  USE lsode_quantities ! B^Phi = B*grad PHi = B^sPhi_s + B^u Phi_u + B^vPhi_v
7  IMPLICIT NONE ! ... and inverts to get B^s and B^u.
8  INTEGER, INTENT(IN):: neq ! Recall that Phi = v, thus Phi_s=Phi_u = 0
9  REAL(rprec), INTENT(IN):: phi
10  REAL(rprec), DIMENSION(neq), INTENT(IN) :: yin
11  REAL(rprec), DIMENSION(neq), INTENT(OUT):: fout
12  INTEGER:: i, isize
13  REAL(KIND=rprec) :: s_eval, phi_eval, theta_eval,
14  & bs, bu, bv, jac
15 
16  isize = neq/2
17  DO i = 1, isize !loop over field lines
18  s_eval = yin(2*i-1); theta_eval = yin(2*i)
19  phi_eval = phi
20  CALL reset_coords(s_eval, theta_eval, phi_eval)
21  s_nb(i) = s_eval; phi_nb(i) = phi_eval
22  theta_nb(i) = theta_eval
23  ENDDO
24 
25  CALL dbs3vd(0, 0, 0, isize, s_nb, theta_nb, phi_nb, ks, ku, ! Br
26  & kv, s_knots, u_knots, v_knots, is, itht, iphi, br_coef,
27  & br_s_nb, is_uniformx, is_uniformy, is_uniformz )
28 
29  CALL dbs3vd(0, 0, 0, isize, s_nb, theta_nb, phi_nb, ks, ku, ! Bphi
30  & kv, s_knots, u_knots, v_knots, is, itht, iphi, bphi_coef,
31  & bphi_s_nb, is_uniformx, is_uniformy, is_uniformz )
32 
33  CALL dbs3vd(0, 0, 0, isize, s_nb, theta_nb, phi_nb, ks, ku, ! Bz
34  & kv, s_knots, u_knots, v_knots, is, itht, iphi, bz_coef,
35  & bz_s_nb, is_uniformx, is_uniformy, is_uniformz )
36 
37  CALL dbs3vd(0, 0, 0, isize, s_nb, theta_nb, phi_nb, ks, ku, ! R
38  & kv, s_knots, u_knots, v_knots, is, itht, iphi, r_coef,
39  & r_s_nb, is_uniformx, is_uniformy, is_uniformz )
40 
41  DO i = 1, isize
42  s_eval = s_nb(i); theta_eval = theta_nb(i)
43  phi_eval = phi_nb(i)
44  drds_s_nb(i) = dbs3dr(1, 0, 0, s_eval, theta_eval, phi_eval, ! dR/ds
45  & ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi, r_coef)
46  drdtht_s_nb(i) = dbs3dr(0, 1, 0, s_eval, theta_eval, phi_eval, ! dR/du
47  & ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi, r_coef)
48  drdphi_s_nb(i) = dbs3dr(0, 0, 1, s_eval, theta_eval, phi_eval, ! dR/dv
49  & ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi, r_coef)
50  dzds_s_nb(i) = dbs3dr(1, 0, 0, s_eval, theta_eval, phi_eval, ! dZ/ds
51  & ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi, z_coef)
52  dzdtht_s_nb(i) = dbs3dr(0, 1, 0, s_eval, theta_eval, phi_eval, ! dZ/du
53  & ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi, z_coef)
54  dzdphi_s_nb(i) = dbs3dr(0, 0, 1, s_eval, theta_eval, phi_eval, ! dZ/dv
55  & ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi, z_coef)
56  ENDDO
57 
58  DO i = 1, isize ! Field line ODEs in cylindrical geometry
59  jac = drds_s_nb(i)*dzdtht_s_nb(i) - drdtht_s_nb(i)*dzds_s_nb(i)
60  bs = r_s_nb(i)*(dzdtht_s_nb(i)*br_s_nb(i) - !B^s = (R(Zu*B^R-RuB^Z))+B^Phi*(RuZv-RvZu))/Jacob
61  & drdtht_s_nb(i)*bz_s_nb(i)) + bphi_s_nb(i)*(
62  & drdtht_s_nb(i)*dzdphi_s_nb(i) - drdphi_s_nb(i)*dzdtht_s_nb(i))
63  bu = r_s_nb(i)*(drds_s_nb(i)*bz_s_nb(i) - !B^u = (R(RsB^z-ZsB^R)}B^Phi*(ZsRv-RsZu))/Jacob
64  & dzds_s_nb(i)*br_s_nb(i)) + bphi_s_nb(i)*(
65  & drdphi_s_nb(i)*dzds_s_nb(i) - drds_s_nb(i)*dzdphi_s_nb(i))
66  bv = bphi_s_nb(i) !B^v = B^Phi
67  IF (bv == zero) stop 'BV becomes zero for some orbit!'
68  fout(2*i-1) = bs/(jac*bv) ! r.h.s. = B^s/B^v
69  fout(2*i) = bu/(jac*bv) ! r.h.s. = B^u/B^v
70  END DO
71 
72  END SUBROUTINE rt_hand_side_dphi