rt_hand_side_dphi.f

      SUBROUTINE rt_hand_side_dphi(neq, phi, yin, fout)                         ! Defines field-line equations for LSODE integrator for type="toroidal"
      USE stel_kinds                                                            ! The right hand sides are B^s/B^v  and B^u/B^v.
      USE stel_constants                                                        ! To get them, one simply uses:
      USE bspline                                                               ! B^R =B*grad R = B^sRs + B^uRu + B^vRv
      USE spline3d_fit_coefs                                                    ! B^z =B*grad Z = B^sZs + B^uZu + B^vZu
      USE lsode_quantities                                                      ! B^Phi = B*grad PHi = B^sPhi_s + B^u Phi_u + B^vPhi_v
      IMPLICIT NONE                                                             ! ... and inverts to get B^s and B^u.
      INTEGER, INTENT(IN):: neq                                                 ! Recall that Phi = v, thus Phi_s=Phi_u = 0
      REAL(rprec), INTENT(IN):: phi
      REAL(rprec), DIMENSION(neq), INTENT(IN) :: yin
      REAL(rprec), DIMENSION(neq), INTENT(OUT):: fout
      INTEGER:: i, isize
      REAL(KIND=rprec) :: s_eval, phi_eval, theta_eval,
     &  bs, bu, bv, jac
 
      isize = neq/2
      DO i = 1, isize                                                           !loop over field lines
        s_eval = yin(2*i-1); theta_eval = yin(2*i)
        phi_eval = phi
        CALL reset_coords(s_eval, theta_eval, phi_eval)
        s_nb(i) = s_eval; phi_nb(i) = phi_eval
        theta_nb(i) = theta_eval
      ENDDO
 
      CALL dbs3vd(0, 0, 0, isize, s_nb, theta_nb, phi_nb, ks, ku,               ! Br
     &  kv, s_knots, u_knots, v_knots, is, itht, iphi, br_coef,
     &  br_s_nb, is_uniformx, is_uniformy, is_uniformz )
 
      CALL dbs3vd(0, 0, 0, isize, s_nb, theta_nb, phi_nb, ks, ku,               ! Bphi
     &  kv, s_knots, u_knots, v_knots, is, itht, iphi, bphi_coef,
     &  bphi_s_nb, is_uniformx, is_uniformy, is_uniformz )
 
      CALL dbs3vd(0, 0, 0, isize, s_nb, theta_nb, phi_nb, ks, ku,               ! Bz
     &   kv, s_knots, u_knots, v_knots, is, itht, iphi, bz_coef,
     &   bz_s_nb, is_uniformx, is_uniformy, is_uniformz )
 
      CALL dbs3vd(0, 0, 0, isize, s_nb, theta_nb, phi_nb, ks, ku,               ! R
     &  kv, s_knots, u_knots, v_knots, is, itht, iphi, r_coef,
     &  r_s_nb, is_uniformx, is_uniformy, is_uniformz )
 
      DO i = 1, isize
        s_eval = s_nb(i); theta_eval = theta_nb(i)
        phi_eval = phi_nb(i)
        drds_s_nb(i) = dbs3dr(1, 0, 0, s_eval, theta_eval, phi_eval,            ! dR/ds
     &    ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi, r_coef)
        drdtht_s_nb(i) = dbs3dr(0, 1, 0, s_eval, theta_eval, phi_eval,          ! dR/du
     &    ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi, r_coef)
        drdphi_s_nb(i) = dbs3dr(0, 0, 1, s_eval, theta_eval, phi_eval,          ! dR/dv
     &    ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi, r_coef)
        dzds_s_nb(i) = dbs3dr(1, 0, 0, s_eval, theta_eval, phi_eval,            ! dZ/ds
     &    ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi, z_coef)
        dzdtht_s_nb(i) = dbs3dr(0, 1, 0, s_eval, theta_eval, phi_eval,          ! dZ/du
     &    ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi, z_coef)
        dzdphi_s_nb(i) = dbs3dr(0, 0, 1, s_eval, theta_eval, phi_eval,          ! dZ/dv
     &    ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi, z_coef)
      ENDDO
 
      DO i = 1, isize                                                           ! Field line ODEs in cylindrical geometry
        jac = drds_s_nb(i)*dzdtht_s_nb(i) - drdtht_s_nb(i)*dzds_s_nb(i)
        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
     &    drdtht_s_nb(i)*bz_s_nb(i)) + bphi_s_nb(i)*(
     &    drdtht_s_nb(i)*dzdphi_s_nb(i) - drdphi_s_nb(i)*dzdtht_s_nb(i))
        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
     &    dzds_s_nb(i)*br_s_nb(i)) + bphi_s_nb(i)*(
     &    drdphi_s_nb(i)*dzds_s_nb(i) - drds_s_nb(i)*dzdphi_s_nb(i))
        bv = bphi_s_nb(i)                                                       !B^v = B^Phi
        IF (bv == zero) stop 'BV becomes zero for some orbit!'
        fout(2*i-1) = bs/(jac*bv)                                               ! r.h.s. = B^s/B^v
        fout(2*i) = bu/(jac*bv)                                                 ! r.h.s. = B^u/B^v
      END DO
 
      END SUBROUTINE rt_hand_side_dphi