diff --git a/src/gray_core.f90 b/src/gray_core.f90
index bb3c2e8..cfc9672 100644
--- a/src/gray_core.f90
+++ b/src/gray_core.f90
@@ -2207,9 +2207,10 @@ bb: do
   subroutine print_bres(bres)
     ! Prints the EC resonance surface table (unit 70)
 
-    use equilibrium,  only : rmnm, rmxm, zmnm, zmxm, bfield
-    use units,        only : ubres, unit_active
-    use magsurf_data, only : npsi
+    use const_and_precisions, only : comp_eps
+    use equilibrium,          only : rmnm, rmxm, zmnm, zmxm, bfield
+    use units,                only : ubres, unit_active
+    use magsurf_data,         only : npsi
 
     ! subroutine arguments
     real(wp_), intent(in) :: bres
@@ -2228,10 +2229,10 @@ bb: do
     if (.not. unit_active(ubres)) return
 
     ! Build a regular (R, z) grid
-    dr = (rmxm - rmnm)/(npsi - 1)
+    dr = (rmxm - rmnm - comp_eps)/(npsi - 1)
     dz = (zmxm - zmnm)/(npsi - 1)
     do j=1,npsi
-      rv(j) = rmnm + dr*(j - 1)
+      rv(j) = comp_eps + rmnm + dr*(j - 1)
       zv(j) = zmnm + dz*(j - 1)
     end do
 
@@ -2271,10 +2272,11 @@ bb: do
   subroutine print_maps(B_res, xgcn, R0, Npl0)
     ! Prints several input quantities on a regular grid (unit 72)
 
-    use equilibrium,  only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield
-    use coreprofiles, only : density, temp
-    use units,        only : umaps, unit_active
-    use magsurf_data, only : npsi
+    use const_and_precisions, only : comp_eps
+    use equilibrium,          only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield
+    use coreprofiles,         only : density, temp
+    use units,                only : umaps, unit_active
+    use magsurf_data,         only : npsi
 
     ! subroutine arguments
     real(wp_), intent(in) :: B_res  ! resonant magnetic field, e/m_eω
@@ -2291,10 +2293,10 @@ bb: do
     if (.not. unit_active(umaps)) return
 
     ! Build a regular (R, z) grid
-    dR = (rmxm - rmnm)/(npsi - 1)
+    dR = (rmxm - rmnm - comp_eps)/(npsi - 1)
     dz = (zmxm - zmnm)/(npsi - 1)
     do j = 1, npsi
-      R(j) = rmnm + dR*(j - 1)
+      R(j) = comp_eps + rmnm + dR*(j - 1)
       z(j) = zmnm + dz*(j - 1)
     end do