diff --git a/src/core_ocean/analysis_members/Registry_eddy_product_variables.xml b/src/core_ocean/analysis_members/Registry_eddy_product_variables.xml
index 48be8b032c..696466c101 100644
--- a/src/core_ocean/analysis_members/Registry_eddy_product_variables.xml
+++ b/src/core_ocean/analysis_members/Registry_eddy_product_variables.xml
@@ -39,12 +39,28 @@
 		<var name="velocityMeridionalTimesTemperature" type="real" dimensions="nVertLevels nCells Time" units="m s^{-1} C"
 			description="cell-wise product of component of horizontal velocity in the northward direction and temperature"
 		/>
-                <var name="velocityZonalTimesTemperature_GM" type="real" dimensions="nVertLevels nCells Time" units="m s^{-1} C"
-                        description="cell-wise product of component of horizontal bolus velocity in the eastward direction and temperature"
-                />
-                <var name="velocityMeridionalTimesTemperature_GM" type="real" dimensions="nVertLevels nCells Time" units="m s^{-1} C"
-                        description="cell-wise product of component of horizontal bolus velocity in the northward direction and temperature"
-                />
+		<var name="normalVelocitySquared" type="real" dimensions="nVertLevels nEdges Time" units="m^2 s^{-2}"
+			description="edge based square of normal velocity"
+		/>
+		<var name="normalVelocityTimesTemperature" type="real" dimensions="nVertLevels nEdges Time" units="m s^{-1} C"
+			description="edge based product of normal velocity and temperature"
+		/>
+		<var name="velocityZonalTimesTemperature_GM" type="real" dimensions="nVertLevels nCells Time" units="m s^{-1} C"
+			description="cell-wise product of component of horizontal bolus velocity in the eastward direction and temperature"
+			 packages="gm"
+		/>
+		<var name="velocityMeridionalTimesTemperature_GM" type="real" dimensions="nVertLevels nCells Time" units="m s^{-1} C"
+			description="cell-wise product of component of horizontal bolus velocity in the northward direction and temperature"
+			 packages="gm"
+		/>
+		<var name="normalGMBolusVelocitySquared" type="real" dimensions="nVertLevels nEdges Time" units="m^2 s^{-2}"
+			description="edge based square of normal velocity"
+			 packages="gm"
+		/>
+		<var name="normalGMBolusVelocityTimesTemperature" type="real" dimensions="nVertLevels nEdges Time" units="m s^{-1} C"
+			description="edge based product of normal velocity and temperature"
+			 packages="gm"
+		/>
 	</var_struct>
 	<streams>
 		<stream name="eddyProductVariablesOutput" type="output"
@@ -62,7 +78,9 @@
 			<var name="velocityMeridionalSquared"/>
 			<var name="velocityZonalTimesTemperature"/>
 			<var name="velocityMeridionalTimesTemperature"/>
-                        <var name="velocityZonalTimesTemperature_GM"/>
-                        <var name="velocityMeridionalTimesTemperature_GM"/>
+			<var name="velocityZonalTimesTemperature_GM"/>
+			<var name="velocityMeridionalTimesTemperature_GM"/>
+			<var name="normalVelocityTimesTemperature"/>
+			<var name="normalVelocitySquared"/>
 		</stream>
 	</streams>
diff --git a/src/core_ocean/analysis_members/mpas_ocn_eddy_product_variables.F b/src/core_ocean/analysis_members/mpas_ocn_eddy_product_variables.F
index a7f7da1666..35f735ebf7 100644
--- a/src/core_ocean/analysis_members/mpas_ocn_eddy_product_variables.F
+++ b/src/core_ocean/analysis_members/mpas_ocn_eddy_product_variables.F
@@ -44,6 +44,7 @@ module ocn_eddy_product_variables
    use mpas_stream_manager
 
    use ocn_constants
+   use ocn_config
    use ocn_diagnostics_routines
 
    implicit none
@@ -177,15 +178,17 @@ subroutine ocn_compute_eddy_product_variables(domain, timeLevel, err)!{{{
       type (mpas_pool_type), pointer :: tracersPool
       type (mpas_pool_type), pointer :: diagnosticsPool
 
-      integer, pointer :: nVertLevels, nCellsSolve
-      integer :: iTracer, k, iCell
+      integer, pointer :: nEdges, nVertLevels, nCellsSolve
+      integer :: iTracer, k, iCell, iEdge, cell1, cell2
       integer, pointer :: index_temperature, index_salinity
-      integer, dimension(:), pointer :: maxLevelCell
+      integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell
+      integer, dimension(:,:), pointer :: cellsOnEdge
 
       real (kind=RKIND), dimension(:), pointer :: ssh, SSHSquared
       real (kind=RKIND), dimension(:,:), pointer :: velocityZonal, velocityMeridional, GMBolusVelocityZonal, GMBolusVelocityMeridional, &
            velocityZonalSquared, velocityMeridionalSquared, velocityZonalTimesTemperature, velocityMeridionalTimesTemperature, &
-           velocityZonalTimesTemperature_GM, velocityMeridionalTimesTemperature_GM
+           velocityZonalTimesTemperature_GM, velocityMeridionalTimesTemperature_GM, normalVelocity, normalVelocitySquared, &
+           normalVelocityTimesTemperature, normalGMBolusVelocity, normalGMBolusVelocityTimesTemperature, normalGMBolusVelocitySquared
       real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers
 
       err = 0
@@ -202,37 +205,107 @@ subroutine ocn_compute_eddy_product_variables(domain, timeLevel, err)!{{{
 
          call mpas_pool_get_dimension(block % dimensions, 'nVertLevels', nVertLevels)
          call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve)
+         call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges)
 
+         call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 1)
          call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 1)
          call mpas_pool_get_dimension(tracersPool, 'index_temperature', index_temperature)
 
          call mpas_pool_get_array(statePool, 'ssh',ssh, 1)
          call mpas_pool_get_array(diagnosticsPool, 'velocityZonal', velocityZonal)
          call mpas_pool_get_array(diagnosticsPool, 'velocityMeridional', velocityMeridional)
-         call mpas_pool_get_array(diagnosticsPool, 'GMBolusVelocityZonal',GMBolusVelocityZonal)
-         call mpas_pool_get_array(diagnosticsPool, 'GMBolusVelocityMeridional', GMBolusVelocityMeridional)
          call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
+         call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
+         call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop)
+
          call mpas_pool_get_array(eddyProductVariablesAMPool, 'SSHSquared', SSHSquared)
          call mpas_pool_get_array(eddyProductVariablesAMPool, 'velocityZonalSquared', velocityZonalSquared)
          call mpas_pool_get_array(eddyProductVariablesAMPool, 'velocityMeridionalSquared', velocityMeridionalSquared)
          call mpas_pool_get_array(eddyProductVariablesAMPool, 'velocityZonalTimesTemperature', velocityZonalTimesTemperature)
          call mpas_pool_get_array(eddyProductVariablesAMPool, 'velocityMeridionalTimesTemperature', &
             velocityMeridionalTimesTemperature)
-         call mpas_pool_get_array(eddyProductVariablesAMPool,'velocityZonalTimesTemperature_GM', velocityZonalTimesTemperature_GM)
-         call mpas_pool_get_array(eddyProductVariablesAMPool,'velocityMeridionalTimesTemperature_GM', &
-            velocityMeridionalTimesTemperature_GM)
-
-         do iCell = 1,nCellsSolve
-            SSHSquared(iCell) = ssh(iCell)**2
-            do k = 1, maxLevelCell(iCell)
-               velocityZonalSquared(k,iCell) = velocityZonal(k,iCell)**2
-               velocityMeridionalSquared(k,iCell) = velocityMeridional(k,iCell)**2
-               velocityZonalTimesTemperature(k,iCell) = velocityZonal(k,iCell)*activeTracers(index_temperature,k,iCell)
-               velocityMeridionalTimesTemperature(k,iCell) = velocityMeridional(k,iCell)*activeTracers(index_temperature,k,iCell)
-               velocityZonalTimesTemperature_GM(k,iCell) = GMBolusVelocityZonal(k,iCell)*activeTracers(index_temperature,k,iCell)
-               velocityMeridionalTimesTemperature_GM(k,iCell) = GMBolusVelocityMeridional(k,iCell)*activeTracers(index_temperature,k,iCell)
-             end do
-         end do
+         call mpas_pool_get_array(eddyProductVariablesAMPool,'normalVelocityTimesTemperature', normalVelocityTimesTemperature)
+         call mpas_pool_get_array(eddyProductVariablesAMPool,'normalVelocitySquared', normalVelocitySquared)
+
+         ! if GM is active, include the GM variables.
+         ! I repeated the block of code here for better performance, rather than
+         ! split the GM and non-GM variables into separate loops.
+         if (config_use_GM) then
+            call mpas_pool_get_array(diagnosticsPool, 'normalGMBolusVelocity', normalGMBolusVelocity)
+            call mpas_pool_get_array(diagnosticsPool, 'GMBolusVelocityZonal',GMBolusVelocityZonal)
+            call mpas_pool_get_array(diagnosticsPool, 'GMBolusVelocityMeridional', GMBolusVelocityMeridional)
+            call mpas_pool_get_array(eddyProductVariablesAMPool,'velocityZonalTimesTemperature_GM', velocityZonalTimesTemperature_GM)
+            call mpas_pool_get_array(eddyProductVariablesAMPool,'velocityMeridionalTimesTemperature_GM', &
+               velocityMeridionalTimesTemperature_GM)
+            call mpas_pool_get_array(eddyProductVariablesAMPool,'normalGMBolusVelocityTimesTemperature', normalGMBolusVelocityTimesTemperature)
+            call mpas_pool_get_array(eddyProductVariablesAMPool,'normalGMBolusVelocitySquared', normalGMBolusVelocitySquared)
+
+            !$omp parallel
+            !$omp do schedule(runtime) private(cell1, cell2, k)
+            do iEdge = 1,nEdges
+               cell1 = cellsOnEdge(1,iEdge)
+               cell2 = cellsOnEdge(2,iEdge)
+
+               do k = 1, maxLevelEdgeTop(iEdge)
+                  normalVelocityTimesTemperature(k,iEdge) = 0.5_RKIND*(activeTracers(index_temperature,k,cell1) + &
+                               activeTracers(index_temperature,k,cell2)) * normalVelocity(k,iEdge)
+                  normalVelocitySquared(k,iEdge) = normalVelocity(k,iEdge)**2
+                  normalGMBolusVelocityTimesTemperature(k,iEdge) = 0.5_RKIND*(activeTracers(index_temperature,k,cell1) + &
+                               activeTracers(index_temperature,k,cell2)) * normalGMBolusVelocity(k,iEdge)
+                  normalGMBolusVelocitySquared(k,iEdge) = normalGMBolusVelocity(k,iEdge)**2
+               end do
+            end do
+            !$omp end do
+            !$omp end parallel
+
+            !$omp parallel
+            !$omp do schedule(runtime) private(k)
+            do iCell = 1,nCellsSolve
+               SSHSquared(iCell) = ssh(iCell)**2
+               do k = 1, maxLevelCell(iCell)
+                  velocityZonalSquared(k,iCell) = velocityZonal(k,iCell)**2
+                  velocityMeridionalSquared(k,iCell) = velocityMeridional(k,iCell)**2
+                  velocityZonalTimesTemperature(k,iCell) = velocityZonal(k,iCell)*activeTracers(index_temperature,k,iCell)
+                  velocityMeridionalTimesTemperature(k,iCell) = velocityMeridional(k,iCell)*activeTracers(index_temperature,k,iCell)
+                  velocityZonalTimesTemperature_GM(k,iCell) = GMBolusVelocityZonal(k,iCell)*activeTracers(index_temperature,k,iCell)
+                  velocityMeridionalTimesTemperature_GM(k,iCell) = GMBolusVelocityMeridional(k,iCell)*activeTracers(index_temperature,k,iCell)
+                end do
+            end do
+            !$omp end do
+            !$omp end parallel
+
+         ! if GM is not active, do exactly the same computation as above but
+         ! without the GM variables.
+         else
+            !$omp parallel
+            !$omp do schedule(runtime) private(cell1, cell2, k)
+            do iEdge = 1,nEdges
+               cell1 = cellsOnEdge(1,iEdge)
+               cell2 = cellsOnEdge(2,iEdge)
+
+               do k = 1, maxLevelEdgeTop(iEdge)
+                  normalVelocityTimesTemperature(k,iEdge) = 0.5_RKIND*(activeTracers(index_temperature,k,cell1) + &
+                               activeTracers(index_temperature,k,cell2)) * normalVelocity(k,iEdge)
+                  normalVelocitySquared(k,iEdge) = normalVelocity(k,iEdge)**2
+               end do
+            end do
+            !$omp end do
+            !$omp end parallel
+
+            !$omp parallel
+            !$omp do schedule(runtime) private(k)
+            do iCell = 1,nCellsSolve
+               SSHSquared(iCell) = ssh(iCell)**2
+               do k = 1, maxLevelCell(iCell)
+                  velocityZonalSquared(k,iCell) = velocityZonal(k,iCell)**2
+                  velocityMeridionalSquared(k,iCell) = velocityMeridional(k,iCell)**2
+                  velocityZonalTimesTemperature(k,iCell) = velocityZonal(k,iCell)*activeTracers(index_temperature,k,iCell)
+                  velocityMeridionalTimesTemperature(k,iCell) = velocityMeridional(k,iCell)*activeTracers(index_temperature,k,iCell)
+                end do
+            end do
+            !$omp end do
+            !$omp end parallel
+         endif
 
          block => block % next
       end do