diff --git a/src/gray_core.f90 b/src/gray_core.f90
index 2106d78..42365b1 100644
--- a/src/gray_core.f90
+++ b/src/gray_core.f90
@@ -142,7 +142,7 @@ contains
     ! print Btot=Bres
     ! print ne, Te, q, Jphi versus psi, rhop, rhot
     call print_bres(bres)
-    call print_prof
+    call print_prof(params%profiles)
     call print_maps(bres, xgcn,                                &
                     norm2(params%antenna%pos(1:2)) * 0.01_wp_, &
                     sin(params%antenna%beta*degree))
@@ -2180,27 +2180,48 @@ bb: do
   end subroutine print_headers
 
 
-  subroutine print_prof
-    ! Prints the (input) plasma kinetic profiles (unit 98)
+  subroutine print_prof(params)
+    ! Prints the (input) plasma kinetic profiles (unit 55)
 
+    use gray_params,   only : profiles_parameters
     use equilibrium,   only : q_spline, fq, frhotor, tor_curr_psi
     use coreprofiles,  only : density, temp
     use units,         only : uprfin, unit_active
 
     implicit none
 
+    ! suborutine arguments
+    type(profiles_parameters), intent(in) :: params
+
     ! local constants
     real(wp_), parameter :: eps = 1.e-4_wp_
 
     ! local variables
-    integer :: i
-    real(wp_) :: psin, rhot, jphi, dens, ddens
+    integer :: i, N_data, N_tail
+    real(wp_) :: rhot, jphi, dens, ddens
+    real(wp_) :: psin, dpsin, psin_last
 
     if (.not. unit_active(uprfin)) return
 
+    ! N of profiles data points
+    N_data = q_spline%ndata
+
+    ! parameters for the density tail (numerical profiles only)
+    if (params%iprof == 1) then
+      N_tail = N_data / 10                          ! N of density tail points
+      psin_last = q_spline%data(N_data)             ! last data point
+      dpsin = (params%psnbnd - psin_last) / N_tail  ! Δψ for uniform grid
+    else
+      N_tail = 0
+    end if
+
     write (uprfin, *) '#psi rhot ne Te q Jphi'
-    do i = 1, q_spline%ndata
-      psin = q_spline%data(i)
+    do i = 1, N_data + N_tail
+      if (i > N_data) then
+        psin = psin_last + dpsin*(i - N_data)
+      else
+        psin = q_spline%data(i)
+      end if
       rhot = frhotor(sqrt(psin))
       call density(psin, dens, ddens)
       jphi = tor_curr_psi(max(eps, psin))
diff --git a/src/main.f90 b/src/main.f90
index c1444d5..20293b8 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -502,7 +502,7 @@ contains
 
     ! Print ne, Te, q, Jphi versus psi, rhop, rhot
     call print_bres(bres)
-    call print_prof
+    call print_prof(params%profiles)
     call print_maps(bres, xgcn,                                &
                     norm2(params%antenna%pos(1:2)) *0.01_wp_ , &
                     sin(params%antenna%beta*degree))