diff --git a/components/mpas-albany-landice/driver/glc_comp_mct.F b/components/mpas-albany-landice/driver/glc_comp_mct.F
index 8bc9fc255d17..ba1043e52e7f 100644
--- a/components/mpas-albany-landice/driver/glc_comp_mct.F
+++ b/components/mpas-albany-landice/driver/glc_comp_mct.F
@@ -1410,7 +1410,7 @@ subroutine glc_import_mct(x2g_g, errorCode)
             n = n + 1
             sfcMassBal(i) = x2g_g % rAttr(index_x2g_Flgl_qice, n)
             floatingBasalMassBal(i) = x2g_g % rAttr(index_x2g_Fogx_qiceli, n)
-            ismip6_2dThermalForcing(i) = x2g_g % rAttr(index_x2g_So_tf300, n)
+            ismip6_2dThermalForcing(i) = x2g_g % rAttr(index_x2g_So_tf2d, n)
 !            surfaceTemperature(i) = x2g_g % rAttr(index_x2g_Sl_tsrf, n)
 !JW            basalOceanHeatflx(i) = x2g_g % rAttr(index_x2g_Fogo_qiceh, n)
 !            basalOceanHeatflx(i) = x2g_g % rAttr(index_x2g_Fogx_qicehi, n)
diff --git a/components/mpas-albany-landice/driver/glc_cpl_indices.F b/components/mpas-albany-landice/driver/glc_cpl_indices.F
index e84ac9174b19..459b534a7e07 100644
--- a/components/mpas-albany-landice/driver/glc_cpl_indices.F
+++ b/components/mpas-albany-landice/driver/glc_cpl_indices.F
@@ -22,7 +22,7 @@ module glc_cpl_indices
   integer, public :: index_x2g_So_htv    = 0 !Ice shelf ocean heat transfer velocity
   integer, public :: index_x2g_So_stv    = 0 !Ice shelf ocean salinity transfer velocity
   integer, public :: index_x2g_So_rhoeff = 0 !Ocean effective pressure
-  integer, public :: index_x2g_So_tf300  = 0 !Ocean thermal forcing at 300m
+  integer, public :: index_x2g_So_tf2d   = 0 !Ocean thermal forcing at predefined critical depth
   integer, public :: index_x2g_Fogx_qiceli = 0 !Subshelf mass flux
   integer, public :: index_x2g_Fogx_qicehi = 0 !Subshelf heat flux for the ice sheet
 
@@ -71,7 +71,7 @@ subroutine glc_cpl_indices_set( )
     index_x2g_Fogx_qiceli =               mct_avect_indexra(x2g,'Fogx_qiceli',perrwith='quiet')
     index_x2g_Fogx_qicehi =               mct_avect_indexra(x2g,'Fogx_qicehi',perrwith='quiet')
     index_x2g_So_rhoeff =                 mct_avect_indexra(x2g,'So_rhoeff',perrwith='quiet')
-    index_x2g_So_tf300  =                 mct_avect_indexra(x2g,'So_tf300',perrwith='quiet')
+    index_x2g_So_tf2d =                   mct_avect_indexra(x2g,'So_tf2d',perrwith='quiet')
 
     !Following block of x2g/g2x vectors are used internally within coupler for subshelf melt flux
     !calculations (and so do not have directly-related export-side arrays)
diff --git a/components/mpas-ocean/driver/mpaso_cpl_indices.F b/components/mpas-ocean/driver/mpaso_cpl_indices.F
index c51aa68bf5b8..c5e84d2509d9 100644
--- a/components/mpas-ocean/driver/mpaso_cpl_indices.F
+++ b/components/mpas-ocean/driver/mpaso_cpl_indices.F
@@ -37,7 +37,7 @@ module mpaso_cpl_indices
   integer :: index_o2x_So_htv !ocean heat-transfer velocity
   integer :: index_o2x_So_stv !ocean salt-transfer velocity
   integer :: index_o2x_So_rhoeff !ocean effective density
-  integer :: index_o2x_So_tf300 !ocean thermal forcing at 300m
+  integer :: index_o2x_So_tf2d !ocean thermal forcing at predefined critical depth
 
 
   ! ocn -> drv  (BGC)
@@ -209,7 +209,7 @@ subroutine mpaso_cpl_indices_set( )
     index_o2x_So_htv        = mct_avect_indexra(o2x,'So_htv')
     index_o2x_So_stv        = mct_avect_indexra(o2x,'So_stv')
     index_o2x_So_rhoeff     = mct_avect_indexra(o2x,'So_rhoeff')
-    index_o2x_So_tf300      = mct_avect_indexra(o2x,'So_tf300')
+    index_o2x_So_tf2d       = mct_avect_indexra(o2x,'So_tf2d')
 
     index_o2x_So_algae1     = mct_avect_indexra(o2x,'So_algae1',perrWith='quiet')
     index_o2x_So_algae2     = mct_avect_indexra(o2x,'So_algae2',perrWith='quiet')
diff --git a/components/mpas-ocean/driver/ocn_comp_mct.F b/components/mpas-ocean/driver/ocn_comp_mct.F
index e237360034ad..415c4582c130 100644
--- a/components/mpas-ocean/driver/ocn_comp_mct.F
+++ b/components/mpas-ocean/driver/ocn_comp_mct.F
@@ -2703,7 +2703,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
                                                avgRemovedIceRunoffFlux, &
                                                avgLandIceHeatFlux, &
                                                avgRemovedIceRunoffHeatFlux, &
-                                               avgThermalForcing300m
+                                               avgThermalForcingAtCritDepth
 
    real (kind=RKIND), dimension(:,:), pointer :: avgTracersSurfaceValue, avgSurfaceVelocity, &
                                                  avgSSHGradient, avgOceanSurfacePhytoC, &
@@ -2765,7 +2765,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
      call mpas_pool_get_array(forcingPool, 'avgSurfaceVelocity', avgSurfaceVelocity)
      call mpas_pool_get_array(forcingPool, 'avgSSHGradient', avgSSHGradient)
      call mpas_pool_get_array(forcingPool, 'avgTotalFreshWaterTemperatureFlux', avgTotalFreshWaterTemperatureFlux)
-     call mpas_pool_get_array(forcingPool, 'avgThermalForcing300m', avgThermalForcing300m)
+     call mpas_pool_get_array(forcingPool, 'avgThermalForcingAtCritDepth', avgThermalForcingAtCritDepth)
 
      if ( frazilIceActive ) then
         call mpas_pool_get_array(forcingPool, 'seaIceEnergy', seaIceEnergy)
@@ -2942,7 +2942,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
           o2x_o % rAttr(index_o2x_So_stv, n) = landIceTracerTransferVelocities(indexSaltTrans,i)
           o2x_o % rAttr(index_o2x_So_rhoeff, n) = 0.0_RKIND
        endif
-       o2x_o % rAttr(index_o2x_So_tf300, n) = avgThermalForcing300m(i)
+       o2x_o % rAttr(index_o2x_So_tf2d, n) = avgThermalForcingAtCritDepth(i)
 
 
        !Fyke: test
diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index a04379e50367..2ce561550d16 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -770,6 +770,10 @@
 					description="If true, solid and liquid runoff from the Antarctic Ice Sheet (below 60S latitude) coming from the coupled is zeroed in the coupler import routines.  To be used with data iceberg fluxes coming from the sea ice model."
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_2d_thermal_forcing_depth" type="real" default_value="300" units="m"
+					description="Depth at which to pass 2d thermal forcing to the coupler for use in the GLC component.  Note that mapping files for this field must be created with a mask to exclude ocean grid cells shallower than this value and thus must be regenerated if this value is changed."
+					possible_values="any non-negative value"
+		/>
 	</nml_record>
 	<nml_record name="shortwaveRadiation" mode="init;forward">
 		<nml_option name="config_sw_absorption_type" type="character" default_value="none"
@@ -3896,8 +3900,8 @@
 			description="The time-averaged effective ocean density within ice shelves based on Archimedes' principle."
 			packages="landIceCouplingPKG"
 		/>
-		<var name="avgThermalForcing300m" type="real" dimensions="nCells Time" units="C"
-			description="The time-averaged thermal forcing at 300 m depth"
+		<var name="avgThermalForcingAtCritDepth" type="real" dimensions="nCells Time" units="C"
+			description="The time-averaged thermal forcing at the predefined critical depth specified by config_2d_thermal_forcing_depth"
 		/>
 		<!-- Input fields for data (prescribed) land-ice fluxes -->
 		<var name="dataLandIceFreshwaterFlux" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F b/components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F
index 34d2bd989f50..db3d100b4e2f 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F
@@ -60,7 +60,7 @@ subroutine ocn_time_average_coupled_init(forcingPool)!{{{
                                                     avgRemovedRiverRunoffFlux, avgRemovedIceRunoffFlux, &
                                                     avgLandIceHeatFlux, avgRemovedIceRunoffHeatFlux
 
-        real (kind=RKIND), dimension(:), pointer :: avgThermalForcing300m
+        real (kind=RKIND), dimension(:), pointer :: avgThermalForcingAtCritDepth
 
         integer :: iCell
         integer, pointer :: nAccumulatedCoupled, nCells
@@ -154,12 +154,12 @@ subroutine ocn_time_average_coupled_init(forcingPool)!{{{
            !$omp end parallel
         end if
 
-        call mpas_pool_get_array(forcingPool, 'avgThermalForcing300m', avgThermalForcing300m)
+        call mpas_pool_get_array(forcingPool, 'avgThermalForcingAtCritDepth', avgThermalForcingAtCritDepth)
 
         !$omp parallel
         !$omp do schedule(runtime)
         do iCell = 1, nCells
-           avgThermalForcing300m(iCell) = 0.0_RKIND
+           avgThermalForcingAtCritDepth(iCell) = 0.0_RKIND
         end do
         !$omp end do
         !$omp end parallel
@@ -272,13 +272,13 @@ subroutine ocn_time_average_coupled_accumulate(statePool, forcingPool, timeLevel
                                                     removedRiverRunoffFlux, avgRemovedRiverRunoffFlux, &
                                                     removedIceRunoffFlux, avgRemovedIceRunoffFlux, &
                                                     avgRemovedIceRunoffHeatFlux, &
-                                                    avgThermalForcing300m
+                                                    avgThermalForcingAtCritDepth
 
         type (mpas_pool_type), pointer :: tracersPool
 
         real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers
         integer, pointer :: indexTemperature, indexSalinity
-        integer :: iLevel0300, iLevel
+        integer :: iLevelCritDepth, iLevel
         real (kind=RKIND) :: freezingTemp
 
         real (kind=RKIND), dimension(:,:,:), pointer :: &
@@ -417,17 +417,19 @@ subroutine ocn_time_average_coupled_accumulate(statePool, forcingPool, timeLevel
 
         end if
 
-        call mpas_pool_get_array(forcingPool, 'avgThermalForcing300m', avgThermalForcing300m)
+        call mpas_pool_get_array(forcingPool, 'avgThermalForcingAtCritDepth', avgThermalForcingAtCritDepth)
         call mpas_pool_get_subpool(statePool, 'tracers', tracersPool)
         call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 2)
         call mpas_pool_get_dimension(tracersPool, 'index_temperature', indexTemperature)
         call mpas_pool_get_dimension(tracersPool, 'index_salinity', indexSalinity)
 
-        ! find find vertical level that is just above the 300 m reference level
-        iLevel0300 = 1
+
+        ! find vertical level that is just above the critical depth reference level
+        ! this does not account for depression due to ice shelf cavities or sea ice
+        iLevelCritDepth = 1
         do iLevel = 1, nVertLevels
-           if(refBottomDepth(iLevel) > 300.0_RKIND) then
-              iLevel0300 = iLevel-1
+           if(refBottomDepth(iLevel) > config_2d_thermal_forcing_depth) then
+              iLevelCritDepth = iLevel-1
               exit
            end if
         end do
@@ -436,13 +438,13 @@ subroutine ocn_time_average_coupled_accumulate(statePool, forcingPool, timeLevel
         ! calculate thermal forcing at identified level for each cell
         do iCell = 1, nCells
            ! ignore cells that are too shallow
-           if (iLevel0300 <= maxLevelCell(iCell)) then
+           if (iLevelCritDepth <= maxLevelCell(iCell)) then
               ! this uses the level shallower than the reference level.  could interpolate instead
               ! note: assuming no LandIce cavity, but we may want to support that
-              freezingTemp = ocn_freezing_temperature(salinity=activeTracers(indexSalinity, iLevel0300, iCell), &
-                 pressure=pressure(iLevel0300, iCell),  inLandIceCavity=.false.)
-              avgThermalForcing300m(iCell) = ( avgThermalForcing300m(iCell) * nAccumulatedCoupled &
-                 + activeTracers(indexTemperature, iLevel0300, iCell) - freezingTemp ) / ( nAccumulatedCoupled + 1)
+              freezingTemp = ocn_freezing_temperature(salinity=activeTracers(indexSalinity, iLevelCritDepth, iCell), &
+                 pressure=pressure(iLevelCritDepth, iCell),  inLandIceCavity=.false.)
+              avgThermalForcingAtCritDepth(iCell) = ( avgThermalForcingAtCritDepth(iCell) * nAccumulatedCoupled &
+                 + activeTracers(indexTemperature, iLevelCritDepth, iCell) - freezingTemp ) / ( nAccumulatedCoupled + 1)
             end if
         end do
         !$omp end do
diff --git a/driver-mct/shr/seq_flds_mod.F90 b/driver-mct/shr/seq_flds_mod.F90
index 6ecdd5598fc7..8d7404b8a854 100644
--- a/driver-mct/shr/seq_flds_mod.F90
+++ b/driver-mct/shr/seq_flds_mod.F90
@@ -2987,12 +2987,12 @@ subroutine seq_flds_set(nmlfile, ID, infodata)
        attname  = 'So_rhoeff'
        call metadata_set(attname, longname, stdname, units)
 
-       name = 'So_tf300'
+       name = 'So_tf2d'
        call seq_flds_add(o2x_states,trim(name))
        call seq_flds_add(x2g_states,trim(name))
        call seq_flds_add(x2g_tf_states_from_ocn,trim(name))
-       longname = 'ocean thermal forcing at 300 m depth'
-       stdname  = 'ocean_thermal_forcing_at_300m'
+       longname = 'ocean thermal forcing at predefined critical depth'
+       stdname  = 'ocean_thermal_forcing_at_critical_depth'
        units    = 'C'
        attname  = name
        call metadata_set(attname, longname, stdname, units)