diff --git a/src/core_ocean/Registry.xml b/src/core_ocean/Registry.xml
index 4e1397a381..92fd636456 100644
--- a/src/core_ocean/Registry.xml
+++ b/src/core_ocean/Registry.xml
@@ -391,6 +391,22 @@
 					description="If true, the standard GM for the tracer advection and mixing is turned on. This includes the redi portion which is captured in hmix."
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_gm_lat_variable_c2" type="logical" default_value=".false." units="NA"
+					description="If true, spatially variable phase speed is used Following Ferrari et al (2010)"
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_gm_kappa_lat_depth_variable" type="logical" default_value=".false." units="NA"
+					description="If true, spatitally and deph varying Bolus Kappa is used, following Danabasoglu et al (2007)"
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_gm_min_stratification_ratio" type="real" default_value="0.1" units="NA"
+					description="minimum value for N2/Nmax2 ratio"
+					possible_values="small positive reals"
+		/>
+		<nml_option name="config_gm_min_phase_speed" type="real" default_value="0.1" units="m s^{-1}"
+					description="minimum allowed phase speed for spatially variable computation"
+					possible_values="small positive reals"
+		/>
 		<nml_option name="config_use_Redi_surface_layer_tapering" type="logical" default_value=".false." units="NA"
 					description="If true, the Redi K33 vertical mixing is limited in and just below the ocean boundary layer."
 					possible_values=".true. or .false."
@@ -2549,19 +2565,6 @@
 			 description="kinetic energy of horizonal velocity on cells"
 			 packages="forwardMode;analysisMode"
 		/>
-		<var name="hEddyFlux" type="real" dimensions="nVertLevels nEdges Time" units=""
-			 description="Eddy flux in Gent-McWilliams eddy parameterization"
-			 packages="forwardMode;analysisMode"
-		/>
-		<var name="hKappa" type="real" dimensions="nVertLevels nEdges Time" units=""
-			 description="kappa parameter for Gent-McWilliams eddy parameterization"
-			 packages="forwardMode;analysisMode"
-		/>
-		<var name="hKappaQ" type="real" dimensions="nVertLevels nEdges Time" units=""
-			 description="kappaQ parameter for Gent-McWilliams eddy parameterization"
-			 packages="forwardMode;analysisMode"
-		/>
-
 		<var name="viscosity" type="real" dimensions="nVertLevels nEdges Time" units="m^2 s^{-1}"
 			 description="horizontal viscosity"
 			 packages="forwardMode;analysisMode"
@@ -2667,11 +2670,18 @@
 			 description="Meridional Component of the gradient of sea surface height at cell centers."
 			 packages="forwardMode;analysisMode"
 		/>
-
 		<var name="normalGMBolusVelocity" type="real" dimensions="nVertLevels nEdges Time" units="m s^{-1}"
 			 description="Bolus velocity in Gent-McWilliams eddy parameterization"
 			 packages="forwardMode;analysisMode"
 		/>
+		<var name="cGMphaseSpeed" type="real" dimensions="nEdges Time" units="m s^{-1}"
+			description="phase speed for the bolus velocity calculation"
+			packages="forwardMode;analysisMode"
+			/>
+		<var name="kappaGM3D" type="real" dimensions="nVertLevels nEdges Time" units="m s^{-1}"
+			description="spatially and depth varying GM kappa"
+			packages="forwardMode;analysisMode"
+		/>
 		<var name="GMBolusVelocityX" type="real" dimensions="nVertLevels nCells Time" units="m s^{-1}"
 			 description="Bolus velocity in Gent-McWilliams eddy parameterization, x-direction"
 			 packages="forwardMode;analysisMode"
@@ -3032,7 +3042,6 @@
 		<!-- Coupling fields associated with mass or salinity fluxes -->
 		<var name="evaporationFlux" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
 			 description="Evaporation flux at cell centers from coupler. Positive into the ocean."
-			 packages="thicknessBulkPKG"
 		/>
 		<var name="seaIceSalinityFlux" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
 			 description="Sea ice salinity flux at cell centers from coupler. Positive into the ocean."
@@ -3072,7 +3081,6 @@
 		/>
 		<var name="rainFlux" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
 			 description="Fresh water flux from rain at cell centers from coupler. Positive into the ocean."
-			 packages="thicknessBulkPKG"
 		/>
 		<var name="snowFlux" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
 			 description="Fresh water flux from snow at cell centers from coupler. Positive into the ocean."
@@ -3822,9 +3830,8 @@
 			/>
 	</var_struct>
 
-<!-- Include init modes registry file -->
-#include "mode_init/Registry.xml"
-
+<!-- Include registry files from subdirectories -->
 #include "tracer_groups/Registry_tracers.xml"
 #include "analysis_members/Registry_analysis_members.xml"
+#include "mode_init/Registry.xml"
 </registry>
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 dfbd38e643..48be8b032c 100644
--- a/src/core_ocean/analysis_members/Registry_eddy_product_variables.xml
+++ b/src/core_ocean/analysis_members/Registry_eddy_product_variables.xml
@@ -39,6 +39,12 @@
 		<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_struct>
 	<streams>
 		<stream name="eddyProductVariablesOutput" type="output"
@@ -56,5 +62,7 @@
 			<var name="velocityMeridionalSquared"/>
 			<var name="velocityZonalTimesTemperature"/>
 			<var name="velocityMeridionalTimesTemperature"/>
+                        <var name="velocityZonalTimesTemperature_GM"/>
+                        <var name="velocityMeridionalTimesTemperature_GM"/>
 		</stream>
 	</streams>
diff --git a/src/core_ocean/analysis_members/Registry_moc_streamfunction.xml b/src/core_ocean/analysis_members/Registry_moc_streamfunction.xml
index bff0cabc8c..5824bbf6a0 100644
--- a/src/core_ocean/analysis_members/Registry_moc_streamfunction.xml
+++ b/src/core_ocean/analysis_members/Registry_moc_streamfunction.xml
@@ -63,13 +63,13 @@
 		<var name="mocStreamvalLatAndDepth" type="real" dimensions="nMocStreamfunctionBinsP1 nVertLevels Time" units="sverdrups"
 			description="The value of the MOC streamfunction for each latitude-bin (first dimension) and depth (second dimension)"
 		/>
-		<var name="mocStreamvalLatAndDepthRegion" type="real" dimensions="nMocStreamfunctionBinsP1 nVertLevels R3 Time" units="sverdrups"
+		<var name="mocStreamvalLatAndDepthRegion" type="real" dimensions="nMocStreamfunctionBinsP1 nVertLevels nRegions Time" units="sverdrups"
 			description="The value of the MOC streamfunction for each latitude-bin (first dimension) and depth (second dimension)"
 		/>
 		<var name="binBoundaryMocStreamfunction" type="real" dimensions="nMocStreamfunctionBinsP1" units="varies"
 			 description="Coordinate of southern edge of meridional heat transport bin, either in latitude or y, for plotting."
 		/>
-		<var name="minMaxLatRegion" type="real" dimensions="TWO R3" units="varies"
+		<var name="minMaxLatRegion" type="real" dimensions="TWO nRegions" units="varies"
 			 description="Coordinates of the southern and northern boundaries of each region (for drawing)."
 		/>
 	</var_struct>
diff --git a/src/core_ocean/analysis_members/Registry_time_series_stats_monthly_mean.xml b/src/core_ocean/analysis_members/Registry_time_series_stats_monthly_mean.xml
index 1f91af44de..35336e7203 100644
--- a/src/core_ocean/analysis_members/Registry_time_series_stats_monthly_mean.xml
+++ b/src/core_ocean/analysis_members/Registry_time_series_stats_monthly_mean.xml
@@ -195,5 +195,9 @@
 			<var name="velocityMeridionalSquared"/>
 			<var name="velocityZonalTimesTemperature"/>
 			<var name="velocityMeridionalTimesTemperature"/>
+                        <var name="velocityZonalTimesTemperature_GM"/>
+                        <var name="velocityMeridionalTimesTemperature_GM"/>
+                        <var name="GMBolusVelocityZonal"/>
+                        <var name="GMBolusVelocityMeridional"/>
 		</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 d0b1ea97e6..a7f7da1666 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
@@ -183,8 +183,9 @@ subroutine ocn_compute_eddy_product_variables(domain, timeLevel, err)!{{{
       integer, dimension(:), pointer :: maxLevelCell
 
       real (kind=RKIND), dimension(:), pointer :: ssh, SSHSquared
-      real (kind=RKIND), dimension(:,:), pointer :: velocityZonal, velocityMeridional, &
-           velocityZonalSquared, velocityMeridionalSquared, velocityZonalTimesTemperature, velocityMeridionalTimesTemperature
+      real (kind=RKIND), dimension(:,:), pointer :: velocityZonal, velocityMeridional, GMBolusVelocityZonal, GMBolusVelocityMeridional, &
+           velocityZonalSquared, velocityMeridionalSquared, velocityZonalTimesTemperature, velocityMeridionalTimesTemperature, &
+           velocityZonalTimesTemperature_GM, velocityMeridionalTimesTemperature_GM
       real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers
 
       err = 0
@@ -208,7 +209,8 @@ subroutine ocn_compute_eddy_product_variables(domain, timeLevel, err)!{{{
          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(eddyProductVariablesAMPool, 'SSHSquared', SSHSquared)
          call mpas_pool_get_array(eddyProductVariablesAMPool, 'velocityZonalSquared', velocityZonalSquared)
@@ -216,6 +218,9 @@ subroutine ocn_compute_eddy_product_variables(domain, timeLevel, err)!{{{
          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
@@ -224,7 +229,9 @@ subroutine ocn_compute_eddy_product_variables(domain, timeLevel, err)!{{{
                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
+               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
 
          block => block % next
diff --git a/src/core_ocean/analysis_members/mpas_ocn_moc_streamfunction.F b/src/core_ocean/analysis_members/mpas_ocn_moc_streamfunction.F
index 95828afd5d..38986c4b91 100644
--- a/src/core_ocean/analysis_members/mpas_ocn_moc_streamfunction.F
+++ b/src/core_ocean/analysis_members/mpas_ocn_moc_streamfunction.F
@@ -146,7 +146,7 @@ subroutine ocn_init_moc_streamfunction(domain, err)!{{{
       !! region dimensions
       call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nRegions', nRegions)
       call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nRegionGroups', nRegionGroups)
-      call mpas_pool_get_dimension(domain % blocklist % dimensions, 'R3', maxRegionsInGroup)
+      call mpas_pool_get_dimension(domain % blocklist % dimensions, 'maxRegionsInGroup', maxRegionsInGroup)
 
       !! region config for moc
       call mpas_pool_get_config(domain % configs, 'config_AM_mocStreamfunction_region_group', &
@@ -159,11 +159,19 @@ subroutine ocn_init_moc_streamfunction(domain, err)!{{{
       call mpas_pool_get_array(regionPool, 'regionGroupNames', regionGroupNames)
 
       !!! region preliminaries
+      regionGroupNumber=0
       do i = 1, nRegionGroups
          if (regionGroupNames(i) .eq. additionalRegion) then
             regionGroupNumber = i
          end if
       end do
+      if (regionGroupNumber.eq.0) then
+         call mpas_log_write &
+            ('ocn_moc_streamfunction AM: config_AM_mocStreamfunction_region_group not found in regionGroupNames')
+         call mpas_log_write &
+            ('ocn_moc_streamfunction AM: setting regionGroupNumber=1')
+         regionGroupNumber = 1
+      end if
 
       regionsInAddGroup = MIN(nRegionsInGroup(regionGroupNumber), maxRegionsInGroup)
       !!!! END REGION STUFF
@@ -379,7 +387,7 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{
       !! region dimensions
       call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nRegions', nRegions)
       call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nRegionGroups', nRegionGroups)
-      call mpas_pool_get_dimension(domain % blocklist % dimensions, 'R3', maxRegionsInGroup)
+      call mpas_pool_get_dimension(domain % blocklist % dimensions, 'maxRegionsInGroup', maxRegionsInGroup)
 
       !! region config for moc
       call mpas_pool_get_config(domain % configs, 'config_AM_mocStreamfunction_region_group', &
@@ -392,12 +400,19 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{
       call mpas_pool_get_array(regionPool, 'regionNames', regionNames)
       call mpas_pool_get_array(regionPool, 'regionGroupNames', regionGroupNames)
 
+      regionGroupNumber=0
       do i = 1, nRegionGroups
          if (regionGroupNames(i) .eq. additionalRegion) then
             regionGroupNumber = i
-            !print *, 'found region with the same name:', regionGroupNames(i)
          end if
       end do
+      if (regionGroupNumber.eq.0) then
+         call mpas_log_write &
+            ('ocn_moc_streamfunction AM: config_AM_mocStreamfunction_region_group not found in regionGroupNames')
+         call mpas_log_write &
+            ('ocn_moc_streamfunction AM: setting regionGroupNumber=1')
+         regionGroupNumber = 1
+      end if
 
       regionsInAddGroup = MIN(nRegionsInGroup(regionGroupNumber), maxRegionsInGroup)
 
@@ -426,11 +441,20 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{
       call mpas_pool_get_array(transectPool, 'nTransectsInGroup', nTransectsInGroup)
       call mpas_pool_get_array(transectPool, 'transectNames', transectNames)
       call mpas_pool_get_array(transectPool, 'transectGroupNames', transectGroupNames)
+
+      transectGroupNumber=0
       do i = 1, nTransectGroups
          if (transectGroupNames(i) .eq. additionalTransect) then
             transectGroupNumber = i
          end if
       end do
+      if (transectGroupNumber.eq.0) then
+         call mpas_log_write &
+            ('ocn_moc_streamfunction AM: config_AM_mocStreamfunction_transect_group not found in transectGroupNames')
+         call mpas_log_write &
+            ('ocn_moc_streamfunction AM: setting transectGroupNumber=1')
+         transectGroupNumber = 1
+      end if
 
       transectsInAddGroup = nTransectsInGroup(transectGroupNumber)
 
@@ -493,7 +517,11 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{
          call mpas_pool_get_array(transectPool,'transectEdgeMaskSigns',transectEdgeMaskSigns)
          call mpas_pool_get_array(transectPool,'transectEdgeMasks',transectEdgeMasks)
          call mpas_pool_get_subpool(block % structs, 'state', statePool)
-         call mpas_pool_get_array(statePool, normalVelocityArrayName, normalVelocity, timeLevel)
+         if (trim(normalVelocityArrayName)=="normalVelocity") then
+            call mpas_pool_get_array(statePool, trim(normalVelocityArrayName), normalVelocity, timeLevel)
+         else
+            call mpas_pool_get_array(diagnosticsPool, trim(normalVelocityArrayName), normalVelocity, timeLevel)
+         end if
          call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
          call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop)
          call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
diff --git a/src/core_ocean/shared/mpas_ocn_gm.F b/src/core_ocean/shared/mpas_ocn_gm.F
index 646208519b..1d0d3f9b9f 100644
--- a/src/core_ocean/shared/mpas_ocn_gm.F
+++ b/src/core_ocean/shared/mpas_ocn_gm.F
@@ -51,7 +51,10 @@ module ocn_gm
    logical, pointer :: config_use_Redi_bottom_layer_tapering
    real (kind=RKIND), pointer :: config_Redi_surface_layer_tapering_extent
    real (kind=RKIND), pointer :: config_Redi_bottom_layer_tapering_depth
-
+   logical, pointer :: config_gm_lat_variable_c2
+   logical, pointer :: config_gm_kappa_lat_depth_variable
+   real (kind=RKIND), pointer :: config_gm_min_stratification_ratio
+   real (kind=RKIND), pointer :: config_gm_min_phase_speed
    real (kind=RKIND), parameter :: epsGM = 1.0e-12_RKIND
 
 !***********************************************************************
@@ -98,18 +101,20 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
       !
       !-----------------------------------------------------------------
 
-      real(kind=RKIND), dimension(:,:), pointer :: density, displacedDensity, zMid, normalGMBolusVelocity, hEddyFlux, &
+      real(kind=RKIND), dimension(:,:), pointer :: density, displacedDensity, zMid, normalGMBolusVelocity, &
          layerThicknessEdge, gradDensityEdge, gradDensityTopOfEdge, gradDensityConstZTopOfEdge, gradZMidEdge, &
          gradZMidTopOfEdge, relativeSlopeTopOfEdge, relativeSlopeTopOfCell, k33, gmStreamFuncTopOfEdge, BruntVaisalaFreqTop, &
          gmStreamFuncTopOfCell, dDensityDzTopOfEdge, dDensityDzTopOfCell, relativeSlopeTapering, relativeSlopeTaperingCell, &
-         areaCellSum
-      real(kind=RKIND), dimension(:), pointer   :: boundaryLayerDepth, gmBolusKappa
+         areaCellSum, kappaGM3D
+
+      real(kind=RKIND), dimension(:), pointer   :: boundaryLayerDepth, gmBolusKappa,  cGMphaseSpeed, bottomDepth
       real(kind=RKIND), dimension(:), pointer   :: areaCell, dcEdge, dvEdge, tridiagA, tridiagB, tridiagC, rightHandSide
       integer, dimension(:), pointer   :: maxLevelEdgeTop, maxLevelCell, nEdgesOnCell
       integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
       integer                          :: i, k, iEdge, cell1, cell2, iCell, N, iter
       real(kind=RKIND)                 :: h1, h2, areaEdge, c, BruntVaisalaFreqTopEdge, rtmp, stmp, maxSlopeK33
-
+      real(kind=RKIND)                 :: bottomAv, sumN2, countN2, maxN, kappaSum, ltSum
+      
       ! Dimensions
       integer :: nCells, nEdges
       integer, pointer :: nVertLevels
@@ -124,6 +129,8 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
       call mpas_pool_get_array(diagnosticsPool, 'displacedDensity', displacedDensity)
       call mpas_pool_get_array(diagnosticsPool, 'zMid', zMid)
 
+      call mpas_pool_get_array(diagnosticsPool, 'cGMphaseSpeed', cGMphaseSpeed)
+      call mpas_pool_get_array(diagnosticsPool, 'kappaGM3D', kappaGM3D)
       call mpas_pool_get_array(diagnosticsPool, 'normalGMBolusVelocity', normalGMBolusVelocity)
       call mpas_pool_get_array(diagnosticsPool, 'relativeSlopeTopOfEdge', relativeSlopeTopOfEdge)
       call mpas_pool_get_array(diagnosticsPool, 'relativeSlopeTopOfCell', relativeSlopeTopOfCell)
@@ -131,7 +138,6 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
       call mpas_pool_get_array(diagnosticsPool, 'relativeSlopeTaperingCell', relativeSlopeTaperingCell)
       call mpas_pool_get_array(diagnosticsPool, 'k33', k33)
       call mpas_pool_get_array(diagnosticsPool, 'layerThicknessEdge', layerThicknessEdge)
-      call mpas_pool_get_array(diagnosticsPool, 'hEddyFlux', hEddyFlux)
       call mpas_pool_get_array(diagnosticsPool, 'zMid', zMid)
       call mpas_pool_get_array(diagnosticsPool, 'BruntVaisalaFreqTop', BruntVaisalaFreqTop)
       call mpas_pool_get_array(diagnosticsPool, 'gmStreamFuncTopOfEdge', gmStreamFuncTopOfEdge)
@@ -141,6 +147,7 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
       if (config_use_Redi_surface_layer_tapering) call mpas_pool_get_array(diagnosticsPool, 'boundaryLayerDepth', &
           boundaryLayerDepth)
 
+      call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth)
       call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop',  maxLevelEdgeTop)
       call mpas_pool_get_array(meshPool, 'maxLevelCell',  maxLevelCell)
       call mpas_pool_get_array(meshPool, 'cellsOnEdge',  cellsOnEdge)
@@ -191,22 +198,29 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
 
       ! Assign a huge value to the scratch variables which may manifest itself when
       ! there is a bug.
-      !$omp do schedule(runtime)
-      do iEdge = 1, nEdges + 1
+      !$omp do schedule(runtime) private(k)
+      do iEdge = 1, nEdges
          do k = 1, nVertLevels
             gradDensityEdge(k, iEdge) = huge(0D0)
+            gradZMidEdge(k, iEdge) = huge(0D0)
+            normalGMBolusVelocity(k, iEdge) = 0.0_RKIND
+         end do
+      end do
+      !$omp end do
+
+      !$omp do schedule(runtime) private(k)
+      do iEdge = 1, nEdges
+         do k = 1, nVertLevels + 1
             gradDensityTopOfEdge(k, iEdge) = huge(0D0)
             dDensityDzTopOfEdge(k, iEdge) = huge(0D0)
-            gradZMidEdge(k, iEdge) = huge(0D0)
             gradZMidTopOfEdge(k, iEdge) = huge(0D0)
             relativeSlopeTopOfEdge(k, iEdge) = 0.0_RKIND
             relativeSlopeTapering(k, iEdge) = 0.0_RKIND
-            normalGMBolusVelocity(k, iEdge) = 0.0_RKIND
          end do
       end do
       !$omp end do
 
-      !$omp do schedule(runtime)
+      !$omp do schedule(runtime) private(k)
       do iCell = 1, nCells + 1
          do k = 1, nVertLevels
             dDensityDzTopOfCell(k,  iCell) = huge(0D0)
@@ -516,7 +530,71 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
       !
       !--------------------------------------------------------------------
 
-      c = config_gravWaveSpeed_trunc**2
+      if (config_gm_lat_variable_c2) then
+         !$omp do schedule(runtime) private(cell1, cell2, sumN2, ltSum, countN2, BruntVaisalaFreqTopEdge)
+         do iEdge = 1, nEdges
+            cell1 = cellsOnEdge(1,iEdge)
+            cell2 = cellsOnEdge(2,iEdge)
+            sumN2 = 0.0
+            ltSum = 0.0
+            countN2 = 0
+            
+            do k=2,maxLevelEdgeTop(iEdge)
+
+               BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
+               BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND)
+                
+               sumN2 = sumN2 + BruntVaisalaFreqTopEdge*layerThicknessEdge(k,iEdge)
+               ltSum = ltSum + layerThicknessEdge(k,iEdge)
+               countN2 = countN2 +1
+
+            enddo
+
+            if(countN2 > 0) cGMphaseSpeed(iEdge) = max(config_gm_min_phase_speed ,sqrt(sumN2/ltSum)*ltSum / 3.141592_RKIND)
+
+         enddo
+         !$omp end do
+
+      else
+         !$omp do schedule(runtime)
+         do iEdge = 1, nEdges
+            cGMphaseSpeed(iEdge) = config_gravWaveSpeed_trunc
+         enddo
+         !$omp end do
+      endif
+
+      !$omp do schedule(runtime)
+      do iEdge=1,nEdges
+         kappaGM3D(:,iEdge) = gmBolusKappa(iEdge)
+      enddo 
+      !$omp end do
+
+      if (config_gm_kappa_lat_depth_variable) then
+
+         !$omp do schedule(runtime) private(cell1, cell2, k, BruntVaisalaFreqTopEdge, maxN)
+         do iEdge = 1,nEdges
+            cell1 = cellsOnEdge(1,iEdge)
+            cell2 = cellsOnEdge(2,iEdge)
+
+            maxN = -1.0_RKIND
+            do k=2,maxLevelEdgeTop(iEdge)
+               BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
+               BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND)
+
+               maxN = max(maxN,BruntVaisalaFreqTopEdge)
+
+            enddo
+
+            do k=2,maxLevelEdgeTop(iEdge)
+               BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
+               BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND)
+
+               kappaGM3D(k,iEdge) = gmBolusKappa(iEdge)*max(config_gm_min_stratification_ratio, &
+                       BruntVaisalaFreqTopEdge / (maxN + 1.0E-10_RKIND))
+            enddo
+         enddo
+         !$omp end do
+      endif
 
       nEdges = nEdgesArray( 3 )
 
@@ -533,34 +611,34 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
             k = 2
             BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
             BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND)
-            tridiagB(k-1) = - 2.0_RKIND * config_gravWaveSpeed_trunc**2 / (layerThicknessEdge(k-1,iEdge) &
+            tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1,iEdge) &
                           * layerThicknessEdge(k,iEdge)) - BruntVaisalaFreqTopEdge
-            tridiagC(k-1) = 2.0_RKIND * config_gravWaveSpeed_trunc**2 / layerThicknessEdge(k, iEdge) &
+            tridiagC(k-1) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k, iEdge) &
                           / (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge))
-            rightHandSide(k-1) = gmBolusKappa(iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge)
+            rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge)
 
             ! Second to next to the last rows
             do k = 3, maxLevelEdgeTop(iEdge)-1
                BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
                BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND)
-               tridiagA(k-2) = 2.0_RKIND * config_gravWaveSpeed_trunc**2 / layerThicknessEdge(k-1, iEdge) &
+               tridiagA(k-2) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k-1, iEdge) &
                              / (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge))
-               tridiagB(k-1) = - 2.0_RKIND * config_gravWaveSpeed_trunc**2 / (layerThicknessEdge(k-1, iEdge) &
+               tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1, iEdge) &
                              * layerThicknessEdge(k, iEdge) ) - BruntVaisalaFreqTopEdge
-               tridiagC(k-1) = 2.0_RKIND * config_gravWaveSpeed_trunc**2 / layerThicknessEdge(k, iEdge) &
+               tridiagC(k-1) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k, iEdge) &
                              / (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge))
-               rightHandSide(k-1) = gmBolusKappa(iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge)
+               rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge)
             end do
 
             ! Last row
             k = maxLevelEdgeTop(iEdge)
             BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
             BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND)
-            tridiagA(k-2) = 2.0_RKIND * config_gravWaveSpeed_trunc**2 / layerThicknessEdge(k-1,iEdge) &
+            tridiagA(k-2) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k-1,iEdge) &
                           / (layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge))
-            tridiagB(k-1) = - 2.0_RKIND * config_gravWaveSpeed_trunc**2 / (layerThicknessEdge(k-1, iEdge) &
+            tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1, iEdge) &
                           * layerThicknessEdge(k, iEdge)) - BruntVaisalaFreqTopEdge
-            rightHandSide(k-1) = gmBolusKappa(iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge)
+            rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge)
 
             ! Total number of rows
             N = maxLevelEdgeTop(iEdge) - 1
@@ -747,7 +825,10 @@ subroutine ocn_gm_init(domain, err)!{{{
       call mpas_pool_get_config(ocnConfigs, 'config_use_Redi_bottom_layer_tapering',config_use_Redi_bottom_layer_tapering)
       call mpas_pool_get_config(ocnConfigs, 'config_Redi_surface_layer_tapering_extent',config_Redi_surface_layer_tapering_extent)
       call mpas_pool_get_config(ocnConfigs, 'config_Redi_bottom_layer_tapering_depth',config_Redi_bottom_layer_tapering_depth)
-
+      call mpas_pool_get_config(ocnConfigs, 'config_gm_lat_variable_c2',config_gm_lat_variable_c2)
+      call mpas_pool_get_config(ocnConfigs, 'config_gm_kappa_lat_depth_variable', config_gm_kappa_lat_depth_variable)
+      call mpas_pool_get_config(ocnConfigs, 'config_gm_min_stratification_ratio', config_gm_min_stratification_ratio)
+      call mpas_pool_get_config(ocnConfigs, 'config_gm_min_phase_speed', config_gm_min_phase_speed)
 
       block => domain % blocklist
       do while (associated(block))
@@ -788,6 +869,7 @@ subroutine ocn_gm_init(domain, err)!{{{
             call mpas_dmpar_finalize(domain % dminfo)
          end if
 
+
          block => block % next
       end do
    end subroutine ocn_gm_init!}}}
diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F
index 661e625a6a..4d78f64fa8 100644
--- a/src/core_ocean/shared/mpas_ocn_tendency.F
+++ b/src/core_ocean/shared/mpas_ocn_tendency.F
@@ -117,7 +117,7 @@ subroutine ocn_tend_thick(tendPool, forcingPool, diagnosticsPool, meshPool)!{{{
 
       real (kind=RKIND), dimension(:), pointer :: surfaceThicknessFlux
       real (kind=RKIND), dimension(:), pointer :: surfaceThicknessFluxRunoff
-      real (kind=RKIND), dimension(:,:), pointer :: layerThickness, layerThicknessEdge, &
+      real (kind=RKIND), dimension(:,:), pointer :: layerThicknessEdge, &
          vertAleTransportTop, tend_layerThickness, normalTransportVelocity, fractionAbsorbed, fractionAbsorbedRunoff
 
       integer, pointer :: nCells
@@ -180,7 +180,7 @@ subroutine ocn_tend_thick(tendPool, forcingPool, diagnosticsPool, meshPool)!{{{
       ! surface flux tendency
       !
 
-      call ocn_thick_surface_flux_tend(meshPool, fractionAbsorbed, fractionAbsorbedRunoff, layerThickness,  &
+      call ocn_thick_surface_flux_tend(meshPool, fractionAbsorbed, fractionAbsorbedRunoff, &
          surfaceThicknessFlux, surfaceThicknessFluxRunoff, tend_layerThickness, err)
 
       !
diff --git a/src/core_ocean/shared/mpas_ocn_thick_surface_flux.F b/src/core_ocean/shared/mpas_ocn_thick_surface_flux.F
index d9caec3625..2e1959b414 100644
--- a/src/core_ocean/shared/mpas_ocn_thick_surface_flux.F
+++ b/src/core_ocean/shared/mpas_ocn_thick_surface_flux.F
@@ -72,7 +72,7 @@ module ocn_thick_surface_flux
 !-----------------------------------------------------------------------
 
    subroutine ocn_thick_surface_flux_tend(meshPool, transmissionCoefficients, transmissionCoefficientsRunoff, &
-      layerThickness, surfaceThicknessFlux, surfaceThicknessFluxRunoff, tend, err)!{{{
+      surfaceThicknessFlux, surfaceThicknessFluxRunoff, tend, err)!{{{
       !-----------------------------------------------------------------
       !
       ! input variables
@@ -86,9 +86,6 @@ subroutine ocn_thick_surface_flux_tend(meshPool, transmissionCoefficients, trans
          transmissionCoefficients,     &!< Input: Coefficients for the transmission of surface fluxes
          transmissionCoefficientsRunoff !< Input: Coefficients for the transmission of surface fluxes due to river runoff
 
-      real (kind=RKIND), dimension(:,:), intent(in) :: &
-         layerThickness   !< Input: Layer thickness
-
       real (kind=RKIND), dimension(:), intent(in) :: &
          surfaceThicknessFlux,       &!< Input: surface flux of thickness
          surfaceThicknessFluxRunoff   !< Input: surface flux of thickness due to river runoff
diff --git a/testing_and_setup/compass/README_ocean.md b/testing_and_setup/compass/README_ocean.md
index 1c5ac6d803..5c7d560f61 100644
--- a/testing_and_setup/compass/README_ocean.md
+++ b/testing_and_setup/compass/README_ocean.md
@@ -6,8 +6,8 @@ To set up and run ocean test cases from COMPASS, you will need a conda
 environment.  First, install Miniconda3 (if miniconda is not already
 installed), then create a new conda environment as follows:
 ``` bash
-conda create -n compass_py3.7 -c conda-forge python=3.7 \
-    geometric_features=0.1.1 mpas_tools=0.0.3 jigsaw=0.9.10 jigsawpy=0.0.2 \
+conda create -n compass_py3.7 -c conda-forge -c xylar python=3.7 \
+    geometric_features=0.1.3 mpas_tools=0.0.3 jigsaw=0.9.11 jigsawpy=0.0.2 \
     metis pyflann scikit-image basemap pyamg ffmpeg
 ```
 Each time you want to work with COMPASS, you will need to run:
diff --git a/testing_and_setup/compass/general.config.ocean b/testing_and_setup/compass/general.config.ocean
index 819c41ed1d..d60f520cb4 100644
--- a/testing_and_setup/compass/general.config.ocean
+++ b/testing_and_setup/compass/general.config.ocean
@@ -9,8 +9,8 @@
 # init namelists in the default_inputs directory after a successful build of
 # the ocean model.
 [namelists]
-forward = FULL_PATH_TO_MPAS_MODEL_REPO/namelist.ocean.forward
-init    = FULL_PATH_TO_MPAS_MODEL_REPO/namelist.ocean.init
+forward = FULL_PATH_TO_MPAS_MODEL_REPO/namelist.ocean
+init    = FULL_PATH_TO_MPAS_MODEL_REPO/namelist.ocean
 
 # The streams section defines paths to template streams files that will be used
 # to generate specific streams files. Typically these will point to the forward and
diff --git a/testing_and_setup/compass/ocean/global_ocean/template_forward.xml b/testing_and_setup/compass/ocean/global_ocean/template_forward.xml
index 0af5703b74..5d06bbbf8e 100644
--- a/testing_and_setup/compass/ocean/global_ocean/template_forward.xml
+++ b/testing_and_setup/compass/ocean/global_ocean/template_forward.xml
@@ -16,5 +16,8 @@
 		<option name="config_use_bulk_wind_stress">.true.</option>
 		<option name="config_use_bulk_thickness_flux">.true.</option>
 		<option name="config_use_activeTracers_surface_restoring">.true.</option>
+		<option name="config_use_standardGM">.true.</option>
+		<option name="config_gm_lat_variable_c2">.true.</option>
+		<option name="config_gm_kappa_lat_depth_variable">.true.</option>
 	</namelist>
 </template>