diff --git a/bbb/bbb.v b/bbb/bbb.v
index 387d2d0e..ec28abd2 100755
--- a/bbb/bbb.v
+++ b/bbb/bbb.v
@@ -461,6 +461,11 @@ isnicore(nispmx)  integer   /1,30*0/      +input #switch for ion-density core B.
 				    #=4, use impur. source terms (impur only)
 				    #=5, set d(ni)/dy=-ni/lynicore at midp &
                                     #    ni constant poloidally
+isfniycbozero   real /0./ +input # Switch for divergence-free fluxes on core boundary
+                    #=0, allows divergence-free fluxes to modify net core flux
+                    #=1, redistributes fluxes due to divergence-free term 
+                    #    without affecting the net core boundary flux
+                    #=-1,assumes no divergence-free fluxes on the core boundary
 isupcore(nispmx) integer /nispmx*0/ +input #=0 sets up=upcore on core bdry
 				    #=1 sets d(up)/dy=0 on the core bdry
 				    #=2 sets d^2(up)/dy^2 = 0
diff --git a/bbb/oderhs.m b/bbb/oderhs.m
index c2e2d5e5..9e7ec16d 100755
--- a/bbb/oderhs.m
+++ b/bbb/oderhs.m
@@ -571,7 +571,7 @@ subroutine pandf (xc, yc, neq, time, yl, yldot)
      .     zeffave, noavex, noavey, tiavey, tgavey, psordisold, 
      .     nucxiold(nigmx), nueliold(nigmx), nuelgold(nigmx), rrfac, visxtmp,
      .     vttn, vttp, neavex, pwrebkgold, pwribkgold, feexflr, feixflr,
-     .     naavex,naavey,nuelmolx,nuelmoly
+     .     naavex,naavey,nuelmolx,nuelmoly,fniycboave
       real fqpo, fqpom, friceo, friceom, upeo, upeom, fricio(100), 
      .     friciom(100), upio(100), upiom(100), uupo(100), uupom(100)
       real nevol, ngvol, kionz, krecz, kcxrz, kionm, krecm, kcxrm, nzbg,
@@ -3347,6 +3347,27 @@ ccc          iysptrx is the last closed flux surface (see S.R. nphygeo)
          enddo
       enddo
 
+c ... Normalize core flux to zero to avoid introducing artifical core source/sink
+      if (isfniycbozero .gt. 0) then 
+          fniycboave = 0
+          do ifld = 1, nfsp
+            do ix = ixpt1(1)+1, ixpt2(1)
+              fniycboave = fniycboave + fniycbo(ix, ifld)
+            end do
+            fniycboave = fniycboave / (ixpt2(1) - ixpt1(1))
+            do ix = ixpt1(1)+1, ixpt2(1)
+              fniycbo(ix, ifld) = fniycbo(ix, ifld) - isfniycbozero*fniycboave
+            end do
+          end do
+      else if (isfniycbozero .lt. 0) then 
+          do ifld = 1, nfsp
+            do ix = ixpt1(1)+1, ixpt2(1)
+              fniycbo(ix, ifld) = 0
+            end do
+          end do
+      end if
+             
+
 
 c----------------------------------------------------------------------c
 c          SCALE SOURCE TERMS FROM MONTE-CARLO-NEUTRALS MODEL