minor tweaks to internal kink equilibrium

n-claes · Oct 6, 2021 · 0fcd85f · 0fcd85f
1 parent e601bc8
commit 0fcd85f
Showing 1 changed file with 10 additions and 7 deletions.
diff --git a/src/equilibria/smod_equil_internal_kink_instability.f08 b/src/equilibria/smod_equil_internal_kink_instability.f08
@@ -24,7 +24,7 @@
   module subroutine internal_kink_eq()
     use mod_equilibrium_params, only: cte_rho0, cte_v03, cte_p0, alpha
 
-    real(dp)      :: r, x
+    real(dp)      :: r, x, a0
     real(dp)      :: J0, J1, J2, DJ0, DJ1
     integer       :: i
 
@@ -33,36 +33,39 @@ module subroutine internal_kink_eq()
     )
     call initialise_grid()
 
+    a0 = x_end
     if (use_defaults) then ! LCOV_EXCL_START
       flow = .true.
       cte_rho0 = 1.0d0
       cte_v03  = 1.0d0
       cte_p0 = 9.0d0
-      alpha = 5.0d0 / x_end
+      alpha = 5.0d0 / a0
 
       k2 = 1.0d0
       k3 = 0.16d0 * alpha
     end if ! LCOV_EXCL_STOP
 
     do i = 1, gauss_gridpts
       r = grid_gauss(i)
-      x = r / x_end
+      x = r / a0
 
       J0   = bessel_jn(0, alpha * x)
       J1   = bessel_jn(1, alpha * x)
       J2   = bessel_jn(2, alpha * x)
       DJ0  = -alpha * J1
       DJ1  = alpha * (0.5d0 * J0 - 0.5d0 * J2)
 
-      rho_field % rho0(i) = cte_rho0 * (1.0d0 - x**2)
-      v_field % v03(i)    = cte_v03 * (1.0d0 - x**2)
+      rho_field % rho0(i) = cte_rho0 * (1.0d0 - x**2 / a0**2)
+      v_field % v03(i)    = cte_v03 * (1.0d0 - x**2 / a0**2)
       B_field % B02(i)    = J1
       B_field % B03(i)    = J0
       B_field % B0(i)     = sqrt((B_field % B02(i))**2 + (B_field % B03(i))**2)
       T_field % T0(i)     = cte_p0 / (rho_field % rho0(i))
 
-      rho_field % d_rho0_dr(i) = -2.0d0 * cte_rho0 * x
-      T_field % d_T0_dr(i)     = 2.0d0 * x * cte_p0 / (cte_rho0 * (1.0d0 - x**2)**2)
+      rho_field % d_rho0_dr(i) = -2.0d0 * cte_rho0 * x / a0
+      T_field % d_T0_dr(i)     = 2.0d0 * x * cte_p0 / ( &
+        a0**2 * cte_rho0 * (1.0d0 - x**2)**2 &
+      )
       v_field % d_v03_dr(i)    = -2.0d0 * cte_v03 * x
       B_field % d_B02_dr(i)    = DJ1
       B_field % d_B03_dr(i)    = DJ0