diff --git a/src/jxbforce.f90 b/src/jxbforce.f90
index 6f3345da3..6f49169b8 100644
--- a/src/jxbforce.f90
+++ b/src/jxbforce.f90
@@ -7,7 +7,8 @@
 !> @param bsupv contravariant component of magnetic field \f$B^\zeta\f$
 !> @param bsubu covariant component of magnetic field \f$B_\theta\f$
 !> @param bsubv covariant component of magnetic field \f$B_\zeta\f$
-!> @param bsubsh covariant component of magnetic field \f$B_s\f$ on half grid from bss()
+!> @param bsubsh covariant component of magnetic field \f$B_s\f$: on the half grid
+!>        from bss() on entry, on the full grid on exit
 !> @param bsubsu tangential derivate of covariant component of magnetic field \f$\partial B_s / \partial \theta\f$ (?)
 !> @param bsubsv tangential derivate of covariant component of magnetic field \f$\partial B_s / \partial \zeta\f$  (?)
 !> @param gsqrt Jacobian \f$\sqrt{g}\f$
@@ -33,7 +34,7 @@ SUBROUTINE jxbforce(bsupu, bsupv, bsubu, bsubv, bsubsh, &
   REAL(rprec), DIMENSION(ns,nznt), INTENT(in) :: bsupv
   REAL(rprec), DIMENSION(ns,nznt,0:1), TARGET, INTENT(inout) :: bsubu
   REAL(rprec), DIMENSION(ns,nznt,0:1), TARGET, INTENT(inout) :: bsubv
-  REAL(rprec), DIMENSION(ns,nznt), INTENT(in)  :: bsubsh
+  REAL(rprec), DIMENSION(ns,nznt), INTENT(inout) :: bsubsh
   REAL(rprec), DIMENSION(ns,nznt,0:1) :: bsubsu
   REAL(rprec), DIMENSION(ns,nznt,0:1) :: bsubsv
   REAL(rprec), DIMENSION(ns,nznt), INTENT(in) :: gsqrt
@@ -518,8 +519,13 @@ SUBROUTINE jxbforce(bsupu, bsupv, bsubu, bsubv, bsubsh, &
 
   ! Compute end point values for bsubs
   bsubs(1,:)  = 2*bsubs(2,:)  - bsubs(3,:)
-  !bsubs(ns,:) = 2*bsubs(ns,:) - bsubs(ns-1,:) ! TODO: from ns, ns-1 to ns ???
-  bsubs(ns,:) = 2*bsubs(ns-1,:) - bsubs(ns-2,:)
+  ! The outermost half point sits between bsubs(ns-1) and the boundary, so
+  ! the boundary value follows from it and one full-mesh point.
+  bsubs(ns,:) = 2*bsubsh(ns,:) - bsubs(ns-1,:)
+
+  ! Return the full-mesh B_s to the caller: wrout transforms this array into
+  ! bsubsmns, which the wout file declares on the full mesh.
+  bsubsh = bsubs
 
 
 
diff --git a/src/read_wout_mod.f b/src/read_wout_mod.f
index 7f647f02f..3c1d8dca0 100644
--- a/src/read_wout_mod.f
+++ b/src/read_wout_mod.f
@@ -727,8 +727,8 @@ SUBROUTINE Compute_Currents(bsubsmnc_i, bsubsmns_i,
 !     Computes current harmonics for currXmn == sqrt(g)*JsupX, X = u,v
 !     [Corrected above "JsubX" to "JsupX", JDH 2010-08-16]
 
-!     NOTE: bsub(s,u,v)mn are on HALF radial grid
-!          (in earlier versions, bsubsmn was on FULL radial grid)
+!     NOTE: bsub(u,v)mn are on the HALF radial grid, bsubsmn on the FULL one,
+!          as in ORNL-Fusion/LIBSTELL read_wout_mod.f90
 
 !
       ohs = (ns_i-1)
@@ -744,8 +744,7 @@ SUBROUTINE Compute_Currents(bsubsmnc_i, bsubsmns_i,
 
       DO js = 2, ns_i-1
          WHERE (MOD(INT(xm_nyq_i),2) .EQ. 1)
-            t1 = 0.5_dp*(shalf(js+1)*bsubsmns_i(:,js+1) +
-     &                   shalf(js)  *bsubsmns_i(:,js)) /sfull(js)
+            t1 = bsubsmns_i(:,js)
             bu0 = bsubumnc_i(:,js  )/shalf(js)
             bu1 = bsubumnc_i(:,js+1)/shalf(js+1)
             t2 = ohs*(bu1-bu0)*sfull(js)+0.25_dp*(bu0+bu1)/sfull(js)
@@ -753,7 +752,7 @@ SUBROUTINE Compute_Currents(bsubsmnc_i, bsubsmns_i,
             bv1 = bsubvmnc_i(:,js+1)/shalf(js+1)
             t3 = ohs*(bv1-bv0)*sfull(js)+0.25_dp*(bv0+bv1)/sfull(js)
          ELSEWHERE
-            t1 = 0.5_dp*(bsubsmns_i(:,js+1)+bsubsmns_i(:,js))
+            t1 = bsubsmns_i(:,js)
             t2 = ohs*(bsubumnc_i(:,js+1)-bsubumnc_i(:,js))
             t3 = ohs*(bsubvmnc_i(:,js+1)-bsubvmnc_i(:,js))
          ENDWHERE
@@ -783,8 +782,7 @@ SUBROUTINE Compute_Currents(bsubsmnc_i, bsubsmns_i,
 
       DO js = 2, ns_i-1
          WHERE (MOD(INT(xm_nyq_i),2) .EQ. 1)
-            t1 = 0.5_dp*(shalf(js+1)*bsubsmnc_i(:,js+1)
-     &          +         shalf(js)  *bsubsmnc_i(:,js)) / sfull(js)
+            t1 = bsubsmnc_i(:,js)
             bu0 = bsubumns_i(:,js  )/shalf(js+1)
             bu1 = bsubumns_i(:,js+1)/shalf(js+1)
             t2 = ohs*(bu1-bu0)*sfull(js) + 0.25_dp*(bu0+bu1)/sfull(js)
@@ -792,7 +790,7 @@ SUBROUTINE Compute_Currents(bsubsmnc_i, bsubsmns_i,
             bv1 = bsubvmns_i(:,js+1)/shalf(js+1)
             t3 = ohs*(bv1-bv0)*sfull(js)+0.25_dp*(bv0+bv1)/sfull(js)
          ELSEWHERE
-            t1 = 0.5_dp*(bsubsmnc_i(:,js+1) + bsubsmnc_i(:,js))
+            t1 = bsubsmnc_i(:,js)
             t2 = ohs*(bsubumns_i(:,js+1)-bsubumns_i(:,js))
             t3 = ohs*(bsubvmns_i(:,js+1)-bsubvmns_i(:,js))
          END WHERE
diff --git a/src/wrout.f90 b/src/wrout.f90
index 67c169f67..7b36dfdd8 100644
--- a/src/wrout.f90
+++ b/src/wrout.f90
@@ -602,10 +602,7 @@ SUBROUTINE wrout(bsq, gsqrt, bsubu, bsubv, bsubs, bsupv, bsupu, rzl_array, gc_ar
   bsubumnc(:,1) = 0
   bsubvmnc(:,1) = 0
 
-  ! NOTE: This assumes that bsubs (from which bsubsmns is computed) is on the full grid.
-  ! However, since crmn_o is passed as bsubs in fileout(),
-  ! HERE, bsubs is ACTUALLY on the HALF-grid !!! WTF ???
-  ! (The full-grid bsubs array is a LOCAL variable in jxbforce().)
+  ! bsubs is on the full grid here: jxbforce puts it there.
   bsubsmns(:,1) = 2*bsubsmns(:,2) - bsubsmns(:,3) ! extrapolation on full grid
 
   bsupumnc(:,1) = 0