src/equilibrium.f90: fix J_phi at ψ>1 in print_prof

Since 24edfdc4 print_prof tabulates the input profiles up to ψ_bnd,
however the torr_curr_psi may fail (producing FPEs) outside ψ=1.
Moreover, torr_curr_psi gives the current on the LFS on the specific
line z=z_maxis, while the flux surface average is usually more
interesting, say when comparing J_φ with Jcd.

This change fixes both issues.

src/equilibrium.f90

@@ -65,7 +65,7 @@ module equilibrium
   private
   public read_eqdsk, read_equil_an           ! Reading data files
   public scale_equil, change_cocos           ! Transforming data
-  public fq, bfield, tor_curr, tor_curr_psi  ! Accessing local quantities
+  public fq, bfield, tor_curr                ! Accessing local quantities
   public pol_flux, pol_curr                  !
   public frhotor, frhopol                    ! Converting between poloidal/toroidal flux
   public set_equil_spline, set_equil_an      ! Initialising internal state
@@ -1278,68 +1278,6 @@ contains
   end function tor_curr
 
 
-  function tor_curr_psi(psi_n) result(J_phi)
-    ! Computes the toroidal current J_φ as a function of ψ
-
-    ! function arguments
-    real(wp_), intent(in) :: psi_n
-    real(wp_)             :: J_phi
-
-    ! local constants
-    real(wp_), parameter :: eps = 1.e-4_wp_
-
-    ! local variables
-    real(wp_) :: R1, R2
-
-    call psi_raxis(max(eps, psi_n), R1, R2)
-    J_phi = tor_curr(R2, zmaxis)
-  end function tor_curr_psi
-
-
-  subroutine psi_raxis(psin,r1,r2)
-    use gray_params,          only : iequil
-    use dierckx,              only : profil, sproota
-    use logger,               only : log_error
-
-    ! subroutine arguments
-    real(wp_) :: psin,r1,r2
-
-    ! local constants
-    integer, parameter :: mest=4
-
-    ! local variables
-    integer :: iopt,ier,m
-    real(wp_) :: zc,val
-    real(wp_), dimension(mest) :: zeroc
-    real(wp_), dimension(psi_spline%nknots_x) :: czc
-    character(64) :: msg
-
-    if (iequil < 2) then
-      ! Analytical model
-      val = sqrt(psin)
-      r1 = model%R0 - val*model%a
-      r2 = model%R0 + val*model%a
-    else
-      ! Numerical data
-      iopt=1
-      zc=zmaxis
-      call profil(iopt, psi_spline%knots_x, psi_spline%nknots_x, &
-                  psi_spline%knots_y, psi_spline%nknots_y,       &
-                  psi_spline%coeffs, kspl, kspl, zc,             &
-                  psi_spline%nknots_x, czc, ier)
-      if (ier > 0) then
-        write (msg, '("profil failed with error ",g0)') ier
-        call log_error(msg, mod='equilibrium', proc='psi_raxis')
-      end if
-
-      call sproota(psin, psi_spline%knots_x, psi_spline%nknots_x, &
-                   czc, zeroc, mest, m, ier)
-      r1=zeroc(1)
-      r2=zeroc(2)
-    end if
-  end subroutine psi_raxis
-
-
   subroutine points_ox(rz,zz,rf,zf,psinvf,info)
     ! Finds the location of the O,X points
     use const_and_precisions, only : comp_eps

src/gray_core.f90

@@ -2174,32 +2174,42 @@ bb: do
     ! Prints the (input) plasma kinetic profiles (unit 55)
 
     use gray_params,   only : profiles_parameters
-    use equilibrium,   only : fq, frhotor, tor_curr_psi
+    use equilibrium,   only : fq, frhotor
     use coreprofiles,  only : density, temp
     use units,         only : uprfin, unit_active
-    use magsurf_data,  only : npsi
+    use magsurf_data,  only : npsi, vajphiav
 
     ! suborutine arguments
     type(profiles_parameters), intent(in) :: params
 
     ! local variables
-    integer :: i
+    integer :: i, ntail
     real(wp_) :: rho_t, J_phi, dens, ddens
-    real(wp_) :: psi_n, dpsi_n
+    real(wp_) :: psi_n, rho_p, drho_p
 
     if (.not. unit_active(uprfin)) return
 
-    ! Δψ for the uniform ψ grid
-    dpsi_n = params%psnbnd / (npsi - 1)
+    ! Δρ_p for the uniform ρ_p grid
+    ! Note: we don't use ψ_n because J_phi has been
+    ! sampled uniformly in ρ_p (see flux_average)
+    drho_p = 1.0_wp_ / (npsi - 1)
+
+    ! extra points to reach ψ=ψ_bnd
+    ntail = ceiling((sqrt(params%psnbnd) - 1) / drho_p)
 
     write (uprfin, *) '#psi rhot ne Te q Jphi'
-    do i = 0, npsi - 1
-      psi_n = dpsi_n * i
-      rho_t = frhotor(sqrt(psi_n))
-      J_phi = tor_curr_psi(psi_n)
+    do i = 0, npsi + ntail
+      rho_p = i * drho_p
+      rho_t = frhotor(rho_p)
+      psi_n = rho_p**2
+      if (psi_n < 1) then
+        J_phi = vajphiav(i+1) * 1.e-6_wp_
+      else
+        J_phi = 0
+      end if
       call density(psi_n, dens, ddens)
       write (uprfin, '(12(1x,e12.5))') &
-        psi_n, rho_t, dens, temp(psi_n), fq(psi_n), J_phi*1.e-6_wp_
+        psi_n, rho_t, dens, temp(psi_n), fq(psi_n), J_phi
     end do
   end subroutine print_prof