diff --git a/bbb/bbb.v b/bbb/bbb.v
index 02fe7074..1a3118cf 100755
--- a/bbb/bbb.v
+++ b/bbb/bbb.v
@@ -1669,6 +1669,8 @@ voljcsor(0:nx+1,0:ny+1) _real [A]      #uniform core-region curr sor. in ix,iy
 volpsorg(0:nx+1,0:ny+1,1:ngsp) _real [1/s]#curr source for gas in cell ix,iy
 pondpot(0:nx+1,0:ny+1)  _real [V] /0./ #elec ponderomotive potential
 psgov_use(0:nx+1,0:ny+1,1:ngsp) _real [1/m**3 s]#user-specified gas source
+pwre_use(0:nx+1,0:ny+1)  _real [W/m^3]  #user-given electron energy source, ###-MVU 18-oct-2023
+pwri_use(0:nx+1,0:ny+1)  _real [W/m^3]  #user-given ion energy source, ###-MVU 18-oct-2023
 jcvsor                  real [A]  /0./ #total core-region current for voljcsor
 ix_sjcsor	        integer   /0/  #if nonzero, beginning ix for voljcsor
 ix_ejcsor	        integer   /0/  #if nonzero, ending ix for voljcsor
diff --git a/bbb/oderhs.m b/bbb/oderhs.m
index b8537046..0b0c34a5 100755
--- a/bbb/oderhs.m
+++ b/bbb/oderhs.m
@@ -4732,6 +4732,15 @@ cc         elseif (ishosor .ne. 0)
 
          if (istimingon .eq. 1) call timimpfj (tsimp, xc)
       endif  #loop for isimpon==2
+
+ccc MVU - 18-oct-2023, adding external heating terms
+      do iy=0,ny+1
+         do ix=0,nx+1
+            resee(ix,iy) = resee(ix,iy) + pwre_use(ix,iy)*vol(ix,iy)
+            resei(ix,iy) = resei(ix,iy) + pwri_use(ix,iy)*vol(ix,iy)
+         enddo
+      enddo
+ccc MVU
   
 *  -- joule heating --