diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist
index 19418582863e..e35c24f8be26 100755
--- a/components/mpas-ocean/bld/build-namelist
+++ b/components/mpas-ocean/bld/build-namelist
@@ -1005,6 +1005,12 @@ add_default($nl, 'config_manufactured_solution_wavelength_x');
 add_default($nl, 'config_manufactured_solution_wavelength_y');
 add_default($nl, 'config_manufactured_solution_amplitude');
 
+#####################################
+# Namelist group: init_mode_subgrid #
+#####################################
+
+add_default($nl, 'config_subgrid_table_levels');
+
 ################################################
 # Namelist group: tracer_forcing_activeTracers #
 ################################################
@@ -1775,6 +1781,7 @@ my @groups = qw(run_modes
                 transport_tests
                 init_mode_vert_levels
                 manufactured_solution
+                init_mode_subgrid
                 tracer_forcing_activetracers
                 tracer_forcing_debugtracers
                 tracer_forcing_ecosystracers
diff --git a/components/mpas-ocean/bld/build-namelist-group-list b/components/mpas-ocean/bld/build-namelist-group-list
index 11864939025d..c109ba9990e3 100644
--- a/components/mpas-ocean/bld/build-namelist-group-list
+++ b/components/mpas-ocean/bld/build-namelist-group-list
@@ -40,6 +40,7 @@ my @groups = qw(run_modes
                 transport_tests
                 init_mode_vert_levels
                 manufactured_solution
+                init_mode_subgrid
                 tracer_forcing_activetracers
                 tracer_forcing_debugtracers
                 tracer_forcing_ecosystracers
diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section
index 0c834f05cec1..7acfc0c82307 100644
--- a/components/mpas-ocean/bld/build-namelist-section
+++ b/components/mpas-ocean/bld/build-namelist-section
@@ -522,6 +522,12 @@ add_default($nl, 'config_manufactured_solution_wavelength_x');
 add_default($nl, 'config_manufactured_solution_wavelength_y');
 add_default($nl, 'config_manufactured_solution_amplitude');
 
+#####################################
+# Namelist group: init_mode_subgrid #
+#####################################
+
+add_default($nl, 'config_subgrid_table_levels');
+
 ################################################
 # Namelist group: tracer_forcing_activeTracers #
 ################################################
diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
index 5a3359cd2d7a..1a72035aba3b 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
@@ -646,6 +646,9 @@
 <config_manufactured_solution_wavelength_y>2000000.0</config_manufactured_solution_wavelength_y>
 <config_manufactured_solution_amplitude>1</config_manufactured_solution_amplitude>
 
+<!-- init_mode_subgrid -->
+<config_subgrid_table_levels>-1</config_subgrid_table_levels>
+
 <!-- tracer_forcing_activeTracers -->
 <config_use_activeTracers>.true.</config_use_activeTracers>
 <config_use_activeTracers_surface_bulk_forcing>.true.</config_use_activeTracers_surface_bulk_forcing>
diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
index e01c2b9716f4..4be3335c7aab 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
@@ -2640,6 +2640,17 @@ Default: Defined in namelist_defaults.xml
 </entry>
 
 
+<!-- init_mode_subgrid -->
+
+<entry id="config_subgrid_table_levels" type="integer"
+	category="init_mode_subgrid" group="init_mode_subgrid">
+Number of levels in subgrid lookup tables
+
+Valid values: Any positive non-zero integer. A value of -1 causes this to be overwritten with the configurations subgrid table levels definition.
+Default: Defined in namelist_defaults.xml
+</entry>
+
+
 <!-- tracer_forcing_activeTracers -->
 
 <entry id="config_use_activeTracers" type="logical"
diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index 905343f3e408..f4481a617dfe 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -110,6 +110,9 @@
 			description="Number of wavenumbers used to reconstruct Stokes drift depth profile"
 			definition="6"
 		/>
+		<dim name="nSubgridTableLevels" definition="namelist:config_subgrid_table_levels"
+			description="Number of subgrid look-up table levels"
+		/>
 	</dims>
 
 	<!--*************************************-->
@@ -1183,6 +1186,14 @@
 					description="Safety factor on minimum cell height to ensure the minimum height is not violated due to round-off."
 					possible_values="Real value greater or equal to 0."
 		/>
+		<nml_option name="config_use_subgrid_wetting_drying" type="logical" default_value=".false." units="unitless"
+					description="If true, use subgrid wetting and drying algorithm"
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_use_ssh_gradient_wetting_drying" type="logical" default_value=".false." units="unitless"
+					description="If true, use take into account the ssh gradient across previously flux-limited edges in the wetting and drying algorithm"
+					possible_values=".true. or .false."
+		/>
 	</nml_record>
 	<nml_record name="ocean_constants">
 		<nml_option name="config_density0" type="real" default_value="1026.0" units="kg m^-3"
@@ -1606,6 +1617,12 @@
 					possible_values="Any positive real number"
 		/>
 	</nml_record>
+	<nml_record name="init_mode_subgrid">
+		<nml_option name="config_subgrid_table_levels" type="integer" default_value="-1"
+					description="Number of levels in subgrid lookup tables"
+					possible_values="Any positive non-zero integer. A value of -1 causes this to be overwritten with the configurations subgrid table levels definition."
+		/>
+	</nml_record>
 	<packages>
 		<package name="timeVaryingAtmosphericForcingPKG" description="This package includes variables required for time varying atmospheric forcing"/>
 		<package name="timeVaryingLandIceForcingPKG" description="This package includes variavles required for time varying land-ice forcing"/>
@@ -1637,6 +1654,7 @@
 		<package name="gotmPKG" description="This package is for GOTM variables, which are only needed when using GOTM for the vertical turbulence closure."/>
 		<package name="verticalRemapPKG" description="This package is for variables required for vertical Lagrangian remapping."/>
 		<package name="activeWavePKG" description="This package controls variables required for wave coupling"/>
+		<package name="subgridWetDryPKG" description="This package includes variables required for subgrid wetting and drying."/>
 	</packages>
 
 	<streams>
@@ -2612,6 +2630,67 @@
 		<var name="distanceToCoast" type="real" dimensions="nCells" units="m"
 			 description="Distance of each cell to nearest coastline cell."
 		/>
+                <!-- FIELDS FOR SUBGIRD WETTING AND DRYING-->
+		<var name="subgridWetVolumeCellTable" type="real" dimensions="nSubgridTableLevels nCells" units="m"
+			 description="Pre-computed wet volume per unit area for each cell at given ssh value"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridWetVolumeEdgeTable" type="real" dimensions="nSubgridTableLevels nEdges" units="m"
+			 description="Pre-computed wet volume per unit area for each edge at given ssh value"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridWetVolumeVertexTable" type="real" dimensions="nSubgridTableLevels nVertices" units="m"
+			 description="Pre-computed wet volume per unit area for each edge at given ssh value"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridSshCellTableRange" type="real" dimensions="R3 nCells" units="m"
+			 description="Range of ssh values used for cell lookup table"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridSshEdgeTableRange" type="real" dimensions="R3 nEdges" units="m"
+			 description="Range of ssh vaules used for edge lookup table"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridSshVertexTableRange" type="real" dimensions="R3 nVertices" units="m"
+			 description="Range of ssh vaules used for vertex lookup table"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridEdgeBathymetryMean" type="real" dimensions="nEdges" units="m"
+			 description="Mean bathymetry on edge from subgrid information"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridVertexBathymetryMean" type="real" dimensions="nVertices" units="m"
+			 description="Mean bathymetry on vertex from subgrid information"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridCellBathymetryMin" type="real" dimensions="nCells" units="m"
+			 description="Minimum bathymetry on cell from subgrid information"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridEdgeBathymetryMin" type="real" dimensions="nEdges" units="m"
+			 description="Minimum bathymetry on edge from subgrid information"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridVertexBathymetryMin" type="real" dimensions="nVertices" units="m"
+			 description="Minimum bathymetry on vertex from subgrid information"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridLayerThicknessDebug" type="real" dimensions="nCells" units="m"
+			 description="Debug field"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridWetFractionCellTable" type="real" dimensions="nSubgridTableLevels nCells" units="m"
+			 description="Pre-computed wet fraction for each cell at given ssh value"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridWetFractionEdgeTable" type="real" dimensions="nSubgridTableLevels nEdges" units="m"
+			 description="Pre-computed wet fraction for each edge at given ssh value"
+			 packages="subgridWetDryPKG"
+		/>
+		<var name="subgridWetFractionVertexTable" type="real" dimensions="nSubgridTableLevels nVertices" units="m"
+			 description="Pre-computed wet fraction for each edge at given ssh value"
+			 packages="subgridWetDryPKG"
+		/>
 	</var_struct>
 	<var_struct name="verticalMesh" time_levs="1">
 		<var name="restingThickness" type="real" dimensions="nVertLevels nCells" units="m"
@@ -3436,6 +3515,14 @@
                         description="Coriolis term at the previous time step used in the split-explicit AB2 time stepping"
 			packages="splitAB2TimeIntegrator"
 		/>
+		<var name="layerThicknessCellWetDry" type="real" dimensions="nVertLevels nCells Time" units="m"
+			description="intermediate estimate total layer thickness for wetting/drying purposes"
+			packages="forwardMode"
+		/>
+		<var name="sshCellWetDry" type="real" dimensions="nCells Time" units="m"
+			description="intermediate estimate ssh for wetting/drying purposes"
+			packages="forwardMode"
+		/>
 	</var_struct>
 	<var_struct name="shortwave" time_levs="1">
 		<!-- **********************************************************************
diff --git a/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F b/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F
index ec6a19c81d11..7d9eb083f907 100644
--- a/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F
+++ b/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F
@@ -138,6 +138,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       logical, pointer :: gotmPKGActive
       logical, pointer :: verticalRemapPKGActive
       logical, pointer :: activeWavePKGActive
+      logical, pointer :: subgridWetDryPKGActive
 
       type (mpas_pool_iterator_type) :: pkgItr
       logical, pointer :: packageActive
@@ -172,6 +173,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       logical, pointer :: config_use_time_varying_land_ice_forcing
       logical, pointer :: config_use_gotm
       logical, pointer :: config_use_active_wave
+      logical, pointer :: config_use_subgrid_wetting_drying
 
       character (len=StrKIND), pointer :: config_time_integrator
       character (len=StrKIND), pointer :: config_ocean_run_mode
@@ -427,6 +429,15 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
          activeWavePKGActive = .true.
       endif
 
+      !
+      ! test for use of subgrid wetting and drying
+      !
+      call mpas_pool_get_package(packagePool, 'subgridWetDryPKGActive', subgridWetDryPKGActive)
+      call mpas_pool_get_config(configPool, 'config_use_subgrid_wetting_drying', config_use_subgrid_wetting_drying)
+      if (config_use_subgrid_wetting_drying) then
+         subgridWetDryPKGActive = .true.
+      end if
+
       !
       ! call into analysis member driver to set analysis member packages
       !
diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F
index daf0417a2091..d5ba6b08eb6b 100644
--- a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F
+++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F
@@ -80,6 +80,7 @@ module ocn_forward_mode
    use ocn_submesoscale_eddies
    use ocn_stokes_drift
    use ocn_manufactured_solution
+   use ocn_subgrid
 
    use ocn_high_freq_thickness_hmix_del2
 
@@ -300,6 +301,8 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{
       ierr = ior(ierr, err_tmp)
       call ocn_diagnostics_init(domain, err_tmp)
       ierr = ior(ierr,err_tmp)
+      call ocn_subgrid_init(domain, err_tmp)
+      ierr = ior(ierr,err_tmp)
 
       if(ierr.eq.1) then
           call mpas_log_write('An error was encountered while initializing the MPAS-Ocean forward mode', MPAS_LOG_CRIT)
diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F
index d4f4021eb67f..5fbb56e6ec1f 100644
--- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F
+++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F
@@ -44,6 +44,9 @@ module ocn_time_integration_rk4
    use ocn_surface_land_ice_fluxes
    use ocn_transport_tests
 
+   use ocn_subgrid
+
+
    implicit none
    private
    save
@@ -195,6 +198,9 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
 
       real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersCur, activeTracersNew
 
+      real(kind=RKIND), dimension(:,:), pointer :: subgridSshCellTableRange, subgridWetVolumeCellTable
+
+
       ! Get config options
       call mpas_pool_get_config(domain % configs, 'config_mom_del4', config_mom_del4)
       call mpas_pool_get_config(domain % configs, 'config_filter_btr_mode', config_filter_btr_mode)
@@ -248,6 +254,13 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
          call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
          call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop)
 
+         call mpas_pool_get_array(meshPool, 'subgridWetVolumeCellTable', &
+                                subgridWetVolumeCellTable)
+
+         call mpas_pool_get_array(meshPool, 'subgridSshCellTableRange', &
+                                subgridSshCellTableRange)
+
+
          ! Lower k-loop limit of 1 rather than minLevel* needed in *New = *Cur
          ! assignments below are needed to maintain bit-for-bit results
 
@@ -449,7 +462,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
              if (config_prevent_drying) then
                block => domain % blocklist
                  do while (associated(block))
-                   call ocn_prevent_drying_rk4(block, dt, rk_substep_weights(rk_step), config_zero_drying_velocity, err)
+                   call ocn_prevent_drying_rk4(domain, block, dt, rk_substep_weights(rk_step), config_zero_drying_velocity, err)
                    block => block % next
                  end do
                ! exchange fields for parallelization
@@ -676,6 +689,15 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
                  layerThicknessNew(k, iCell) = tidalInputMask(iCell)*(tidalBCValue(iCell) + bottomDepth(iCell))*(restingThickness(k,iCell)/totalDepth)
                  !(1.0_RKIND - tidalInputMask(iCell))*layerThicknessNew(k, iCell)  ! generalized tappered assumption code
                end do
+
+               if ( config_use_subgrid_wetting_drying ) then
+                  call ocn_subgrid_layer_thickness_lookup( tidalBCValue(iCell), & 
+                                         subgridWetVolumeCellTable(:,iCell), &
+                                         subgridSshCellTableRange(:,iCell),&
+                                         bottomDepth(iCell), &
+                                         layerThicknessNew(1,iCell) )
+               end if
+
              end if
            end do
          end if
diff --git a/components/mpas-ocean/src/mode_init/Makefile b/components/mpas-ocean/src/mode_init/Makefile
index d1aba81e14f5..2f16f1973183 100644
--- a/components/mpas-ocean/src/mode_init/Makefile
+++ b/components/mpas-ocean/src/mode_init/Makefile
@@ -7,7 +7,8 @@ UTILS = mpas_ocn_init_spherical_utils.o \
         mpas_ocn_init_cell_markers.o \
         mpas_ocn_init_interpolation.o \
         mpas_ocn_init_ssh_and_landIcePressure.o \
-        mpas_ocn_init_smoothing.o
+        mpas_ocn_init_smoothing.o \
+        mpas_ocn_init_subgrid.o 
 
 TEST_CASES = mpas_ocn_init_baroclinic_channel.o \
              mpas_ocn_init_lock_exchange.o \
@@ -31,7 +32,8 @@ TEST_CASES = mpas_ocn_init_baroclinic_channel.o \
              mpas_ocn_init_mixed_layer_eddy.o \
              mpas_ocn_init_transport_tests.o \
              mpas_ocn_init_test_sht.o \
-             mpas_ocn_init_parabolic_bowl.o
+             mpas_ocn_init_parabolic_bowl.o \
+	     mpas_ocn_init_buttermilk_bay.o
              #mpas_ocn_init_TEMPLATE.o
 
 all: init_mode
@@ -54,6 +56,8 @@ mpas_ocn_init_vertical_grids.o:
 
 mpas_ocn_init_seaSurfaceHeightAndPressure.o:
 
+mpas_ocn_init_subgrid.o:
+
 mpas_ocn_init_baroclinic_channel.o: $(UTILS)
 
 mpas_ocn_init_iso.o: $(UTILS)
@@ -100,6 +104,7 @@ mpas_ocn_init_test_sht.o: $(UTILS)
 
 mpas_ocn_init_parabolic_bowl.o: $(UTILS)
 
+mpas_ocn_init_buttermilk_bay.o: $(UTILS)
 
 #mpas_ocn_init_TEMPLATE.o: $(UTILS)
 
diff --git a/components/mpas-ocean/src/mode_init/Registry.xml b/components/mpas-ocean/src/mode_init/Registry.xml
index 3c05635a9c1b..1ef457fe5328 100644
--- a/components/mpas-ocean/src/mode_init/Registry.xml
+++ b/components/mpas-ocean/src/mode_init/Registry.xml
@@ -20,6 +20,7 @@
 #include "Registry_mixed_layer_eddy.xml"
 #include "Registry_test_sht.xml"
 #include "Registry_parabolic_bowl.xml"
+#include "Registry_buttermilk_bay.xml"
 // #include "Registry_TEMPLATE.xml"
 
 
diff --git a/components/mpas-ocean/src/mode_init/Registry_buttermilk_bay.xml b/components/mpas-ocean/src/mode_init/Registry_buttermilk_bay.xml
new file mode 100644
index 000000000000..09754d2a7d4b
--- /dev/null
+++ b/components/mpas-ocean/src/mode_init/Registry_buttermilk_bay.xml
@@ -0,0 +1,58 @@
+	<nml_record name="buttermilk_bay" mode="init" configuration="buttermilk_bay">
+		<nml_option name="config_Buttermilk_bay_vert_levels" type="integer" default_value="1" units="unitless"
+			description="Number of vertical levels in Buttermilk bay problem."
+			possible_values="Any positive integer number greater than 1."
+	        />
+	        <nml_option name="config_Buttermilk_bay_adjust_domain" type="logical" default_value=".true." units="unitless"
+			description="Flag to adjust mesh coordinates"
+			possible_values='.true. or .false'
+		/>
+		<nml_option name="config_Buttermilk_bay_subgrid_refinement_level" type="integer" default_value="3" units="unitless"
+			description="Level of refinement to evaluate subgrid bathymety"
+			possible_values="Any positive integer number"
+	        />
+		<nml_option name="config_Buttermilk_bay_subgrid_edge_bathymetry_max_pixel" type="logical" default_value=".true." units="unitless"
+			description="Use maximum elevation pixels as subgrid bathymetry on cell edges"
+			possible_values=".true. or .false."
+	        />
+		<nml_option name="config_Buttermilk_bay_subgrid_table_levels" type="integer" default_value="10" units="unitless"
+			description="Number of levels in suvgrid lookup tables"
+			possible_values="Any positive integer number"
+		/>
+		<nml_option name="config_Buttermilk_bay_subgrid_use_thin_layer" type="logical" default_value=".false." units="unitless"
+			description="Whether to include thin layer in subgrid tables"
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_Buttermilk_bay_topography_source" type="character" default_value="xy_file" units="unitless"
+			description="If 'latlon_file/xy_file', reads in topography from file specified in config_Buttermilk_bay_topography_file"
+			possible_values="'latlon_file' or 'xy_file'"
+		/>
+		<nml_option name="config_Buttermilk_bay_topography_file" type="character" default_value="none" units="unitless"
+			description="Path to the topography initial condition file."
+			possible_values="path/to/topography/file.nc"
+		/>
+		<nml_option name="config_Buttermilk_bay_topography_latlon_degrees"  type="logical" default_value=".false." units="unitless"
+			description="Logical flag that controls if the Lat/Lon fields for topography should be converted to radians. True means input is degrees, false means input is radians."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_Buttermilk_bay_topography_nlat_dimname" type="character" default_value="y" units="unitless"
+			description="Dimension name for the latitude in the topography file."
+			possible_values="Dimension name from input file."
+		/>
+		<nml_option name="config_Buttermilk_bay_topography_nlon_dimname" type="character" default_value="x" units="unitless"
+			description="Dimension name for the longitude in the topography file."
+			possible_values="Dimension name from input file."
+		/>
+		<nml_option name="config_Buttermilk_bay_topography_lat_varname" type="character" default_value="y" units="unitless"
+			description="Variable name for the latitude in the topography file."
+			possible_values="Variable name from input file."
+		/>
+		<nml_option name="config_Buttermilk_bay_topography_lon_varname" type="character" default_value="x" units="unitless"
+			description="Variable name for the longitude in the topography file."
+			possible_values="Variable name from input file."
+		/>
+		<nml_option name="config_Buttermilk_bay_topography_varname" type="character" default_value="Band1" units="unitless"
+			description="Variable name for the topography in the topography file."
+			possible_values="Variable name from input file."
+		/>
+	</nml_record>
diff --git a/components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml b/components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml
index b2188231dc4e..bd01e21e1be7 100644
--- a/components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml
+++ b/components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml
@@ -2,16 +2,16 @@
 		<nml_option name="config_parabolic_bowl_vert_levels" type="integer" default_value="3" units="unitless"
 			description="Number of vertical levels in parabolic bowl."
 			possible_values="Any positive integer"
-	        />
-	        <nml_option name="config_parabolic_bowl_Coriolis_parameter" type="real" default_value="1.031e-4" units="1/s"
+		/>
+		<nml_option name="config_parabolic_bowl_Coriolis_parameter" type="real" default_value="1.031e-4" units="1/s"
 			description="Coriolis paramter"
 			possible_values="Any real number"
 		/>
-	        <nml_option name="config_parabolic_bowl_eta0" type="real" default_value="2.0" units="m"
+		<nml_option name="config_parabolic_bowl_eta0" type="real" default_value="2.0" units="m"
 			description="surface elevation magnitude"
 			possible_values="Any real number"
 		/>
-                <nml_option name="config_parabolic_bowl_b0" type="real" default_value="50.0" units="m"
+		<nml_option name="config_parabolic_bowl_b0" type="real" default_value="50.0" units="m"
 			description="maximum water depth"
 			possible_values="Any real number"
 		/>
@@ -23,13 +23,24 @@
 			description="Gravitational accerlation"
 			possible_values="Any real number"
 		/>
-	        <nml_option name="config_parabolic_bowl_adjust_domain_center" type="logical" default_value="true" units="unitless"
+		<nml_option name="config_parabolic_bowl_adjust_domain_center" type="logical" default_value=".true." units="unitless"
 			description="Flag to adjust mesh coordinates"
+			possible_values='.true. or .false'
+		/>
+		<nml_option name="config_parabolic_bowl_subgrid_refinement_level" type="integer" default_value="3" units="unitless"
+			description="Level of refinement to evaluate subgrid bathymety"
+			possible_values="Any positive integer number"
+		/>
+		<nml_option name="config_parabolic_bowl_subgrid_edge_bathymetry_max_pixel" type="logical" default_value=".true." units="unitless"
+			description="Whether edge subgrid calculations are based on edge CV average or subgrid pixels along edge"
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_parabolic_bowl_subgrid_table_levels" type="integer" default_value="10" units="unitless"
+			description="Number of levels in suvgrid lookup tables"
+			possible_values="Any positive integer number"
+		/>
+		<nml_option name="config_parabolic_bowl_subgrid_use_thin_layer" type="logical" default_value=".false." units="unitless"
+			description="Whether to include thin layer in subgrid tables"
 			possible_values=".true. or .false."
 		/>
 	</nml_record>
-
-  
-
-
-
diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F
new file mode 100644
index 000000000000..0b9f4cbf512b
--- /dev/null
+++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_buttermilk_bay.F
@@ -0,0 +1,906 @@
+! Copyright (c) 2013,  Los Alamos National Security, LLC (LANS)
+! and the University Corporation for Atmospheric Research (UCAR).
+!
+! Unless noted otherwise source code is licensed under the BSD license.
+! Additional copyright and license information can be found in the LICENSE file
+! distributed with this code, or at http://mpas-dev.github.com/license.html
+!
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  ocn_init_buttermilk_bay
+!
+!> \brief  initialize the Buttermilk Bay case
+!> \author D. Wirasaet, S. Brus
+!> \date   May-June 2022
+!> \details
+!>  This module contains the routines for initializing the Buttermilk
+!>  Bay test case with or without subgrid corrections
+!>
+!-----------------------------------------------------------------------
+
+module ocn_init_Buttermilk_bay
+
+   use mpas_kind_types
+   use mpas_io_units
+   use mpas_derived_types
+   use mpas_pool_routines
+   use mpas_constants
+   use mpas_dmpar
+
+   use ocn_constants
+   use ocn_config
+   use ocn_init_vertical_grids
+   use ocn_init_cell_markers
+   use ocn_subgrid
+   use ocn_init_subgrid
+
+   use mpas_constants
+   use mpas_io
+   use mpas_io_streams
+   use mpas_stream_manager
+
+
+   implicit none
+   private
+   save
+
+   !--------------------------------------------------------------------
+   !
+   ! Public parameters
+   !
+   !--------------------------------------------------------------------
+
+   !--------------------------------------------------------------------
+   !
+   ! Public member functions
+   !
+   !--------------------------------------------------------------------
+
+   public :: ocn_init_setup_Buttermilk_bay, &
+             ocn_init_validate_Buttermilk_bay
+
+   !--------------------------------------------------------------------
+   !
+   ! Private module variables
+   !
+   !--------------------------------------------------------------------
+
+   ! For netcdf topobathy input variables
+   integer :: nLatTopo, nLonTopo
+
+   type (field1DReal) :: topoLat, topoLon
+   type (field2DReal) :: topoIC
+
+   real(kind=RKIND), parameter:: eps = 1.0e-10_RKIND ;
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+!  routine ocn_init_setup_Buttermilk_bay
+!
+!> \brief   Setup for this initial condition
+!> \author  D. Wirasaet and S. Brus
+!> \date    May-June 2022
+!> \details
+!>  This routine sets up the initial conditions for this case.
+!>   To be run in sigma vertical coordinates and single-layer
+!
+!-----------------------------------------------------------------------
+
+  subroutine ocn_init_setup_Buttermilk_bay(domain, iErr)!{{{
+    use mpas_vector_operations  ! To calcutate edgeNormalVector
+
+    implicit none
+    !--------------------------------------------------------------------
+
+    type (domain_type), intent(inout) :: domain
+    integer, intent(out) :: iErr
+
+    type (block_type), pointer :: block_ptr
+    type (mpas_pool_type), pointer :: meshPool
+    type (mpas_pool_type), pointer :: statePool
+    type (mpas_pool_type), pointer :: tracersPool
+    type (mpas_pool_type), pointer :: verticalMeshPool
+    type (mpas_pool_type), pointer :: forcingPool
+
+    ! local variables
+    integer :: iCell, iEdge, iVertex, k, idx
+    real (kind=RKIND) :: yMin, yMax, xMin, xMax, dcEdgeMin, dcEdgeMinGlobal
+    real (kind=RKIND) :: yMinGlobal, yMaxGlobal, yMidGlobal, xMinGlobal, xMaxGlobal
+    real (kind=RKIND) :: localVar1, localVar2
+    real (kind=RKIND), dimension(:), pointer :: interfaceLocations
+
+    ! Define dimension pointers
+    integer, pointer :: nCellsSolve, nEdgesSolve, nVerticesSolve, nVertLevels, nVertLevelsP1
+    integer, pointer :: index_temperature, index_salinity
+    integer, pointer :: maxEdges
+
+    ! Define variable pointers
+    logical, pointer :: on_a_sphere
+    integer, dimension(:), pointer :: minLevelCell, maxLevelCell
+    integer, dimension(:), pointer :: nEdgesOnCell
+    integer, dimension(:,:), pointer :: verticesOnCell, verticesOnEdge
+    integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex
+    real (kind=RKIND), dimension(:), pointer :: xCell, yCell, refBottomDepth, refZMid, &
+         vertCoordMovementWeights, bottomDepth, fCell, fEdge, fVertex, dcEdge
+    real (kind=RKIND), dimension(:,:), pointer:: zMid
+
+    real (kind=RKIND), dimension(:), pointer:: xEdge, yEdge, xVertex, yVertex
+    real (kind=RKIND) :: minBottomDepth, maxBottomDepth, globalMaxBottomDepth, globalMinBottomDepth
+    real (kind=RKIND), dimension(:,:), pointer :: layerThickness, restingThickness
+    real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers
+
+    real (kind=RKIND), dimension(:), pointer :: ssh
+    real (kind=RKIND), dimension(:), pointer :: areaCell
+    real (kind=RKIND), dimension(:,:), pointer :: edgeNormalVectors
+    real (kind=RKIND), dimension(:,:), pointer :: normalVelocity
+    ! Elevation Bcs
+    real (kind=RKIND), dimension(:), pointer :: tidalInputMask
+
+
+    real (kind=RKIND):: HH, uu, vv
+    real (kind=RKIND):: RR, num, den
+    real (kind=RKIND):: xshift = 0.0, yshift = 0.0
+    real (kind=RKIND) :: layerThicknessEdgeAverage
+    real (kind=RKIND), dimension(:,:), allocatable :: rSubgridPoints, sSubgridPoints
+    real (kind=RKIND), dimension(:), allocatable :: subgridBathymetryValues, subgridAreas
+    real (kind=RKIND), dimension(:), allocatable :: subgridSshValues
+    real (kind=RKIND), dimension(:), allocatable :: subgridUValues, subgridVValues
+    real (kind=RKIND), dimension(:), allocatable :: uVelocityAverage, vVelocityAverage
+    integer :: nSubgridCell, nSubgridEdge, nSubgridVertex
+    integer :: nSubgridTriPerSlice
+    integer :: v1, v2
+    integer :: c1, c2
+    real (kind=RKIND) :: x(3), y(3)
+    integer :: slice, nSlice
+    real (kind=RKIND) :: deltaZ
+
+
+    integer:: i, j, jj
+    integer:: nsubgridCellEdge, iEdgeSegment
+    real (kind=RKIND):: pi
+    real (kind=RKIND), dimension(:,:), allocatable :: cellEdgeBathymetryValues
+    real (kind=RKIND), dimension(:), allocatable:: dsEdge
+    real (kind=RKIND), dimension(:), allocatable:: xSubgridCell, ySubgridCell
+    real (kind=RKIND):: bathymetryMin, bathymetryMax
+    iErr = 0
+
+    if(config_init_configuration .ne. trim('buttermilk_bay')) return
+
+    ! Determine vertical grid for configuration
+    call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool)
+    call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
+    call mpas_pool_get_dimension(meshPool, 'nVertLevelsP1', nVertLevelsP1)
+    call mpas_pool_get_config(meshPool, 'on_a_sphere', on_a_sphere)
+
+    nVertLevels  = config_Buttermilk_bay_vert_levels ;
+    nVertLevelsP1 = nVertLevels + 1
+
+    allocate(interfaceLocations(nVertLevelsP1))
+    call ocn_generate_vertical_grid( config_vertical_grid, interfaceLocations, ocnConfigs ) ;
+    !! Mental note: interfaceLocatons = (k-1)/N ;
+
+    ! Initalize min/max values to large positive and negative values
+    yMin = 1.0E10_RKIND
+    yMax = -1.0E10_RKIND
+    xMin = 1.0E10_RKIND
+    xMax = -1.0E10_RKIND
+    dcEdgeMin = 1.0E10_RKIND
+
+    ! Determine local min and max values.
+    block_ptr => domain % blocklist
+    do while(associated(block_ptr))
+       call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool)
+
+       call mpas_pool_get_dimension( meshPool, 'nCellsSolve', nCellsSolve )
+       call mpas_pool_get_dimension( meshPool, 'nEdgesSolve', nEdgesSolve )
+
+       call mpas_pool_get_array(meshPool, 'xCell', xCell)
+       call mpas_pool_get_array(meshPool, 'yCell', yCell)
+       call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge)
+
+       yMin = min( yMin, minval(yCell(1:nCellsSolve)))
+       yMax = max( yMax, maxval(yCell(1:nCellsSolve)))
+       xMin = min( xMin, minval(xCell(1:nCellsSolve)))
+       xMax = max( xMax, maxval(xCell(1:nCellsSolve)))
+       dcEdgeMin = min( dcEdgeMin, minval(dcEdge(1:nEdgesSolve)))
+
+       block_ptr => block_ptr % next
+    end do
+
+    ! Determine global min and max values.
+    call mpas_dmpar_min_real(domain % dminfo, yMin, yMinGlobal)
+    call mpas_dmpar_max_real(domain % dminfo, yMax, yMaxGlobal)
+    call mpas_dmpar_min_real(domain % dminfo, xMin, xMinGlobal)
+    call mpas_dmpar_max_real(domain % dminfo, xMax, xMaxGlobal)
+    call mpas_dmpar_min_real(domain % dminfo, dcEdgeMin, dcEdgeMinGlobal)
+
+    pi = acos(-1.0_RKIND)
+
+    xshift = xMin
+    yshift = (3.0_RKIND*yMin + dcEdgeMin*sin(pi/3.0_RKIND))/3.0_RKIND
+    ! print*, "xMin, yMin = ", xMin, yMin, dcEdgeMin
+    ! print*, "xshift, yshift = ", xshift, yshift
+
+    !***********************************************************************
+    !
+    !  Topography
+    !
+    !***********************************************************************
+
+    call mpas_log_write( 'Reading bathymetry from a NetCDF file')
+
+    if (config_Buttermilk_bay_topography_source == 'latlon_file' .or. &
+             config_Buttermilk_bay_topography_source == 'xy_file' ) then
+       call mpas_log_write( 'Reading topography data from file.')
+       call ocn_init_setup_Buttermilk_bay_read_topo(domain, iErr)
+    endif
+
+    !--------------------------------------------------------------------
+    ! Use this section to set initial values
+    !--------------------------------------------------------------------
+
+    block_ptr => domain % blocklist
+    call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool)
+    call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool)
+    call mpas_pool_get_subpool(block_ptr % structs, 'verticalMesh', verticalMeshPool)
+    call mpas_pool_get_subpool(statePool, 'tracers', tracersPool)
+
+    call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool)
+    call mpas_pool_get_array(forcingPool, 'tidalInputMask', tidalInputMask)
+
+    call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) ;
+    call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) ;
+    call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve) ;
+    call mpas_pool_get_dimension(meshPool, 'nVerticesSolve', nVerticesSolve) ;
+    call mpas_pool_get_dimension(meshPool, 'maxEdges', maxEdges)
+
+    call mpas_pool_get_dimension(tracersPool, 'index_temperature', index_temperature)
+    call mpas_pool_get_dimension(tracersPool, 'index_salinity', index_salinity)
+
+    call mpas_pool_get_array(meshPool, 'xCell', xCell)
+    call mpas_pool_get_array(meshPool, 'yCell', yCell)
+    call mpas_pool_get_array(meshPool, 'refBottomDepth', refBottomDepth)
+    call mpas_pool_get_array(meshPool, 'vertCoordMovementWeights', vertCoordMovementWeights)
+    call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth)
+    call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell)
+    call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
+    call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
+    call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
+
+    call mpas_pool_get_array(meshPool, 'verticesOnCell', verticesOnCell)
+    call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge)
+    call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
+    call mpas_pool_get_array(meshPool, 'cellsOnVertex', cellsOnVertex)
+
+    call mpas_pool_get_array(meshPool, 'fCell', fCell)
+    call mpas_pool_get_array(meshPool, 'fEdge', fEdge)
+    call mpas_pool_get_array(meshPool, 'fVertex', fVertex)
+
+    call mpas_pool_get_array(meshPool, 'xEdge', xEdge )
+    call mpas_pool_get_array(meshPool, 'yEdge', yEdge )
+    call mpas_pool_get_array(meshPool, 'xVertex', xVertex )
+    call mpas_pool_get_array(meshPool, 'yVertex', yVertex )
+
+    call mpas_pool_get_array(statePool, 'zMid', zMid, 1) ;
+
+    call mpas_pool_get_array(statePool, 'ssh', ssh, 1)
+    call mpas_pool_get_array(meshPool, 'edgeNormalVectors', edgeNormalVectors ) ;
+    call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity ) ;
+
+    call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 1)
+    call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1)
+
+    call mpas_pool_get_array(verticalMeshPool, 'refZMid', refZMid)
+    call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness)
+
+    call mpas_pool_get_array(meshPool, 'subgridWetVolumeCellTable', &
+                             subgridWetVolumeCellTable)
+    call mpas_pool_get_array(meshPool, 'subgridWetVolumeEdgeTable', &
+                             subgridWetVolumeEdgeTable)
+    call mpas_pool_get_array(meshPool, 'subgridWetVolumeVertexTable', &
+                             subgridWetVolumeVertexTable)
+    call mpas_pool_get_array(meshPool, 'subgridWetFractionCellTable', &
+                             subgridWetFractionCellTable)
+    call mpas_pool_get_array(meshPool, 'subgridWetFractionEdgeTable', &
+                             subgridWetFractionEdgeTable)
+    call mpas_pool_get_array(meshPool,'subgridWetFractionVertexTable',&
+                             subgridWetFractionVertexTable)
+    call mpas_pool_get_array(meshPool, 'subgridSshCellTableRange', &
+                             subgridSshCellTableRange)
+    call mpas_pool_get_array(meshPool, 'subgridSshEdgeTableRange', &
+                             subgridSshEdgeTableRange)
+    call mpas_pool_get_array(meshPool, 'subgridSshVertexTableRange', &
+                             subgridSshVertexTableRange)
+    call mpas_pool_get_array(meshPool, 'subgridEdgeBathymetryMean', &
+                             subgridEdgeBathymetryMean)
+    call mpas_pool_get_array(meshPool, 'subgridVertexBathymetryMean', &
+                             subgridVertexBathymetryMean)
+    call mpas_pool_get_array(meshPool, 'subgridCellBathymetryMin', &
+                             subgridCellBathymetryMin)
+    call mpas_pool_get_array(meshPool, 'subgridEdgeBathymetryMin', &
+                             subgridEdgeBathymetryMin)
+    call mpas_pool_get_array(meshPool, 'subgridVertexBathymetryMin', &
+                             subgridVertexBathymetryMin)
+    call mpas_pool_get_dimension(meshPool, 'nSubgridTableLevels', &
+                                 nSubgridTableLevels)
+    call mpas_pool_get_array(meshPool, 'subgridLayerThicknessDebug', &
+                             subgridLayerThicknessDebug)
+
+    ! if config_buttermilk_bay_adjust_domain_center == .true.,
+    ! Adjust center of the mesh so that its center is located at (0,0)
+    if ( config_Buttermilk_bay_adjust_domain ) then
+       xCell = xCell - xshift ;
+       yCell = yCell - yshift ;
+
+       xEdge = xEdge - xshift ;
+       yEdge = yEdge - yshift ;
+
+       xVertex = xVertex - xshift ;
+       yVertex = yVertex - yshift ;
+
+       ! get min, max coordinates of model domain   !
+       ! after adjusting the coordinates             !
+       yMin = min( yMin, minval(yCell(1:nCellsSolve)))
+       yMax = max( yMax, maxval(yCell(1:nCellsSolve)))
+       xMin = min( xMin, minval(xCell(1:nCellsSolve)))
+       xMax = max( xMax, maxval(xCell(1:nCellsSolve)))
+
+       ! Determine global min and max values.
+       call mpas_dmpar_min_real(domain % dminfo, yMin, yMinGlobal)
+       call mpas_dmpar_max_real(domain % dminfo, yMax, yMaxGlobal)
+       call mpas_dmpar_min_real(domain % dminfo, xMin, xMinGlobal)
+       call mpas_dmpar_max_real(domain % dminfo, xMax, xMaxGlobal)
+    end if
+
+    ! Initlialze vector
+    call mpas_initialize_vectors(meshPool) ;
+
+    minLevelCell(:) = 1
+    do iCell = 1, nCellsSolve
+      ! Set up vertical grid
+      maxLevelCell(iCell) = nVertLevels ; ! sigma coordinates
+    end do
+
+
+    do iCell = 1, nCellsSolve
+       ! Set temperature
+       activeTracers(index_temperature, :, iCell) = 10.0_RKIND
+
+       ! Set salinity
+       activeTracers(index_salinity, :, iCell) = 30.0_RKIND
+
+       ! Set Coriolis parameters, if other than zero
+       fCell(iCell) = 0.0_RKIND;
+    end do
+
+    do iEdge = 1, nEdgesSolve
+       fEdge(iEdge) = 0.0_RKIND;
+    end do
+
+    do iVertex = 1, nVerticesSolve
+       fVertex(iVertex) = 0.0_RKIND;
+    end do
+
+    allocate(uVelocityAverage(nEdgesSolve))
+    allocate(vVelocityAverage(nEdgesSolve))
+
+    if (config_use_subgrid_wetting_drying) then
+
+       call ocn_subgrid_init(domain,iErr)
+       call ocn_init_subgrid_calculations(domain, &
+                                          ocn_init_Buttermilk_bay_bathymetry, &
+                                          ocn_init_Buttermilk_bay_velocity, &
+                                          ocn_init_Buttermilk_bay_ssh, &
+                                          config_Buttermilk_bay_subgrid_refinement_level, &
+                                          config_Buttermilk_bay_subgrid_edge_bathymetry_max_pixel, &
+                                          config_Buttermilk_bay_subgrid_use_thin_layer, &
+                                          uVelocityAverage, &
+                                          vVelocityAverage, &
+                                          iErr)
+    end if
+
+     ! Find max bottom depth
+     maxBottomDepth = maxval( bottomDepth ) ;
+     minBottomDepth = minval( bottomDepth ) ;
+     call mpas_dmpar_max_real( domain % dminfo, maxBottomDepth, globalMaxBottomDepth ) ;
+     call mpas_dmpar_min_real( domain % dminfo, minBottomDepth, globalMinBottomDepth ) ;
+
+     ! Set refBottomDepth and refZMid
+     do k = 1, nVertLevels
+       refBottomDepth(k) = globalMaxBottomDepth*interfaceLocations(k+1) ;
+       refZMid(k) = -0.5_RKIND*( interfaceLocations(k+1) + interfaceLocations(k))*globalMaxBottomDepth ;
+     end do
+
+     ! Set vertCoordMovementWeights
+     vertCoordMovementWeights(:) = 1.0_RKIND
+
+     ! Set velocity
+     do iEdge = 1, nEdgesSolve
+
+       if (config_use_subgrid_wetting_drying) then
+          do k = 1, nVertLevels
+             normalVelocity(k,iEdge) = uVelocityAverage(iEdge)*edgeNormalVectors(1,iEdge) &
+                                     + vVelocityAverage(iEdge)*edgeNormalVectors(2,iEdge) ;
+          end do
+       else
+          call ocn_init_Buttermilk_bay_velocity(xEdge(iEdge), yEdge(iEdge), uu, vv)
+          do k = 1, nVertLevels
+             normalVelocity(k,iEdge) = uu*edgeNormalVectors(1,iEdge) + vv*edgeNormalVectors(2,iEdge) ;
+          end do
+       end if
+     end do
+
+     ! Set layer thickness and ssh
+     if (config_use_wetting_drying) then
+
+       do iCell = 1, nCellsSolve
+         ! Set up vertical grid
+         maxLevelCell(iCell) = nVertLevels ; ! sigma coordinates
+        end do
+
+        do iCell = 1, nCellsSolve
+          !
+          ! make sure depth is thick enough via ssh = TOTAL_DEPTH - bottomDepth
+          ! add a thin layer of nlayer*config_drying_min_cellhight
+          !
+
+          if (config_use_subgrid_wetting_drying) then
+            ! Initial contion for a subgrid run
+            call ocn_subgrid_ssh_lookup(layerThickness(1,iCell),&
+                                        subgridWetVolumeCellTable(:,iCell),&
+                                        subgridSshCellTableRange(:,iCell),&
+                                        bottomDepth(iCell),&
+                                        subgridCellBathymetryMin(iCell),&
+                                        ssh(iCell))
+            !call ocn_subgrid_layer_thickness_lookup(ssh(iCell), &
+            !                               subgridWetVolumeCellTable(:,iCell), &
+            !                               subgridSshCellTableRange(:,iCell),&
+            !                               bottomDepth(iCell),&
+            !                               LayerThickness(1,iCell))
+
+          else
+            ! Initial condition for a standard run
+            call ocn_init_Buttermilk_bay_bathymetry(xCell(iCell), yCell(iCell), bottomDepth(iCell))
+
+            ssh(iCell) = -bottomDepth(iCell) ;
+            do k = 1, maxLevelCell(iCell)
+              !
+              layerThickness(k,iCell) = max( config_drying_min_cell_height, &
+                  bottomDepth(iCell)/real(maxLevelCell(iCell),RKIND) )
+
+
+              if (layerThickness(k,iCell) < config_drying_min_cell_height) then
+                call mpas_log_write('layerThickness($i,$i)=$r', MPAS_LOG_CRIT, &
+                  intArgs=(/k,iCell/), &
+                  realArgs=(/layerThickness(k,iCell)/))
+              end if
+
+              ssh(iCell) = ssh(iCell) + layerThickness(k,iCell) ;
+            end do
+            !
+          endif
+
+          do k = 1, maxLevelCell(iCell)
+            restingThickness(k,iCell) = bottomDepth(iCell)/maxLevelCell(iCell)
+          end do
+
+        end do
+
+     end if
+
+
+     if (config_use_subgrid_wetting_drying) then
+        do iCell = 1,nCellsSolve
+            call ocn_subgrid_layer_thickness_lookup(ssh(iCell), &
+                                                 subgridWetVolumeCellTable(:,iCell), &
+                                                 subgridSshCellTableRange(:,iCell),&
+                                                 bottomDepth(iCell),&
+                                                 subgridLayerThicknessDebug(iCell))
+        enddo
+     end if
+
+     ! Set tidal boundary mask
+     do iCell = 1, nCellsSolve
+       tidalInputMask(iCell) = 0.0_RKIND
+       if ( yCell(iCell) < (yMin+(dcEdgeMin*sin(pi/3.0_RKIND)/2.0_RKIND))  &
+         .and.  yCell(iCell) > (yMin-(dcEDgeMin*sin(pi/3.0_RKIND)/2.0_RKIND)) ) then
+
+          if ( (xCell(iCell) - dcEdgeMin/2.0_RKIND) > 2048.0_RKIND  &
+           .and.  (xCell(iCell) + dcEdgeMin/2.0_RKIND) < 3072.0_RKIND ) then
+              tidalInputMask(iCell) = 1.0_RKIND
+          end if
+        ! spread it over multiple cells
+        ! if (yCell(iCell) > (25.0e3 - 3*dcEdgeMinGlobal)) then
+        !  tidalInputMask(iCell) = exp(-((yCell(iCell)-25.0e3)/dcEdgeMinGlobal)**2.0)
+       end if
+     end do
+
+    deallocate(interfaceLocations)
+    if (config_global_ocean_topography_source == 'latlon_file') then
+        call mpas_log_write( 'Cleaning up topography IC fields')
+        call ocn_init_Buttermilk_bay_destroy_topo_fields()
+    endif
+    !--------------------------------------------------------------------
+
+    print*, "****** End Butter milk bay init *****" ;
+
+    return ;
+  end subroutine ocn_init_setup_Buttermilk_bay!}}}
+
+!***********************************************************************
+!
+!  routine ocn_init_Buttermilk_bay_bathymetry
+!
+!> \brief   Interpolate bathymetry
+!> \author  Steven Brus, D. Wirasaet
+!> \date    November 2022
+!> \details Return the value of the bathymetry at a given x,y point
+!>
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_init_Buttermilk_bay_bathymetry(x, y, depth)!{{{
+
+      implicit none
+
+      real (kind=RKIND), intent(in) :: x, y
+      real (kind=RKIND), intent(out) :: depth
+
+      integer:: ix(2)
+      real (kind=RKIND) :: xc(2)
+      real (kind=RKIND) :: dx(2), x0(2), xN(2), bath(4), f(4), den
+      real (kind=RKIND) :: xv(2), yv(2)
+
+      !--------------------------------------------------------------------
+      !
+      ! DEM must be on a uniform grid
+      x0 = (/ topoLon%array(1), topoLat%array(1) /) ;
+      xN = (/ topoLon%array(nLonTopo), topoLat%array(nLatTopo) /) ;
+      dx = (/ topoLon%array(2) - topoLon%array(1), &
+              topoLat%array(2) - topoLat%array(1) /) ;
+
+      ! Bilienar interpolation
+      xc = (/ x, y /) ;
+      ix = floor( (xc - x0)/dx ) + 1 ;
+
+      if ( ( ix(1) >= 1 .and. ix(1) < nLonTopo ) .and. &
+              ( ix(2) >= 1 .and. ix(2) < nLatTopo ) ) then
+        ! include the west and soutth bourdaries of
+        ! a given uniform raster DEM
+        bath(1) = TopoIC%array(ix(1),ix(2))  ;
+        bath(2) = TopoIC%array(ix(1)+1,ix(2)) ;
+        bath(3) = TopoIC%array(ix(1)+1,ix(2)+1) ;
+        bath(4) = TopoIC%array(ix(1),ix(2)) ;
+
+        xv(1) = topoLon%array(ix(1)) ;
+        xv(2) = topoLon%array(ix(1)+1) ;
+
+        yv(1) = topoLat%array(ix(2)) ;
+        yv(2) = topoLat%array(ix(2)+1) ;
+
+        den = dx(1)*dx(2) ;
+
+        f(1) = ( xc(1) - xv(2) )*( xc(2) - yv(2) )/den ;
+        f(2) = ( xc(1) - xv(1) )*( xc(2) - yv(2) )/(-den) ;
+        f(3) = ( xc(1) - xv(1) )*( xc(2) - yv(1) )/den  ;
+        f(4) = ( xc(1) - xv(2) )*( xc(2) - yv(1) )/(-den) ;
+
+        depth = sum( bath*f ) ;
+      else
+        ! nearest extrapolation !
+        ix(1) = merge( 1, ix(1), ix(1) < 1 ) ;
+        ix(1) = merge( nLonTopo, ix(1), ix(1) >= nLonTopo ) ;
+        ix(2) = merge( 1, ix(2), ix(2) < 1 )  ;
+        ix(2) = merge( nLatTopo, ix(2), ix(2) >= nLatTopo ) ;
+
+        depth = TopoIC%array(ix(1),ix(2)) ;
+      endif
+
+      return ;
+   end subroutine ocn_init_Buttermilk_bay_bathymetry!}}}
+
+!***********************************************************************
+!
+!  routine ocn_init_Buttermilk_bay_ssh
+!
+!> \brief   Compute initial ssh field
+!> \author  Steven Brus, D. Wirasaet
+!> \date    November 2022
+!> \details Use exact solution to compute ssh field for initial conditions
+!>
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_init_Buttermilk_bay_ssh(x, y, bottomDepth, ssh)!{{{
+
+      implicit none
+
+      real (kind=RKIND), intent(in) :: x, y
+      real (kind=RKIND), intent(in) :: bottomDepth
+      real (kind=RKIND), intent(out) :: ssh
+      real (kind=RKIND) :: RR
+
+   !--------------------------------------------------------------------
+
+      ssh = 0.0_RKIND ;
+      ssh = - bottomDepth + max( ssh + bottomDepth, 0.0_RKIND ) ;
+
+      return
+   !--------------------------------------------------------------------
+
+   end subroutine ocn_init_Buttermilk_bay_ssh!}}}
+
+!***********************************************************************
+!
+!  routine ocn_init_Buttermilk_bay_velocity
+!
+!> \brief   Compute initial velocity field
+!> \author  Steven Brus, D. Wirasaet
+!> \date    November 2022
+!> \details Use exact solution to comupte velocity field for initial conditions
+!>
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_init_Buttermilk_bay_velocity(x, y, u, v)!{{{
+
+      implicit none
+
+      real (kind=RKIND), intent(in) :: x, y
+      real (kind=RKIND), intent(out) :: u, v
+      real (kind=RKIND) :: RR, HH
+
+   !--------------------------------------------------------------------
+
+      u = 0.0_RKIND ;
+      v = 0.0_RKIND ;
+   !--------------------------------------------------------------------
+
+   end subroutine ocn_init_Buttermilk_bay_velocity!}}}
+
+!***********************************************************************
+!
+!  routine ocn_init_setup_global_ocean_read_topo
+!
+!> \brief   Read the topography IC file
+!> \author  D. Wirasaet, S. Brus
+!> \date    Sep 2023
+!> \details
+!>  This routine reads the topography IC file, including latitude and longitude
+!>   information for topography data.
+!>
+!>  Adapted from ocn_init_setup_global_ocean_read_topo
+!-----------------------------------------------------------------------
+    subroutine ocn_init_setup_Buttermilk_bay_read_topo(domain, iErr)!{{{
+       type (domain_type), intent(inout) :: domain
+       integer, intent(out) :: iErr
+
+       type (MPAS_Stream_type) :: topographyStream
+
+       iErr = 0
+
+       ! Define stream for depth levels
+       call MPAS_createStream(topographyStream, domain % iocontext, &
+                config_Buttermilk_bay_topography_file, MPAS_IO_NETCDF, &
+                                                   MPAS_IO_READ, ierr=iErr)
+
+       ! Setup topoLat, topoLon, and topoIC fields for stream to be read in
+       topoLat % fieldName = trim(config_Buttermilk_bay_topography_lat_varname)
+       topoLat % dimSizes(1) = nLatTopo
+       topoLat % dimNames(1) = trim(config_Buttermilk_bay_topography_nlat_dimname)
+       topoLat % isVarArray = .false.
+       topoLat % isPersistent = .true.
+       topoLat % isActive = .true.
+       topoLat % hasTimeDimension = .false.
+       topoLat % block => domain % blocklist
+       allocate(topoLat % attLists(1))
+       allocate(topoLat % array(nLatTopo))
+
+       topoLon % fieldName = trim(config_Buttermilk_bay_topography_lon_varname)
+       topoLon % dimSizes(1) = nLonTopo
+       topoLon % dimNames(1) = trim(config_Buttermilk_bay_topography_nlon_dimname)
+       topoLon % isVarArray = .false.
+       topoLon % isPersistent = .true.
+       topoLon % isActive = .true.
+       topoLon % hasTimeDimension = .false.
+       topoLon % block => domain % blocklist
+       allocate(topoLon % attLists(1))
+       allocate(topoLon % array(nLonTopo))
+
+       topoIC % fieldName = trim(config_Buttermilk_bay_topography_varname)
+       topoIC % dimSizes(1) = nLonTopo
+       topoIC % dimSizes(2) = nLatTopo
+       topoIC % dimNames(1) = trim(config_Buttermilk_bay_topography_nlon_dimname)
+       topoIC % dimNames(2) = trim(config_Buttermilk_bay_topography_nlat_dimname)
+       topoIC % isVarArray = .false.
+       topoIC % isPersistent = .true.
+       topoIC % isActive = .true.
+       topoIC % hasTimeDimension = .false.
+       topoIC % block => domain % blocklist
+       allocate(topoIC % attLists(1))
+       allocate(topoIC % array(nLonTopo, nLatTopo))
+
+
+       ! Add topoLat, topoLon, and topoIC fields to stream
+       call MPAS_streamAddField(topographyStream, topoLat, iErr)
+       call MPAS_streamAddField(topographyStream, topoLon, iErr)
+       call MPAS_streamAddField(topographyStream, topoIC, iErr)
+
+       ! Read stream
+       call MPAS_readStream(topographyStream, 1, iErr)
+       topoIC%array = -topoIC%array ;
+
+
+       ! Close stream
+       call MPAS_closeStream(topographyStream)
+
+       if ( config_Buttermilk_bay_topography_latlon_degrees .and. &
+              config_Buttermilk_bay_topography_source == 'latlon_file' ) then
+          topoLat % array(:) = topoLat % array(:) * pii / 180.0_RKIND
+          topoLon % array(:) = topoLon % array(:) * pii / 180.0_RKIND
+       end if
+
+    end subroutine ocn_init_setup_Buttermilk_bay_read_topo!}}}
+
+
+!***********************************************************************
+!
+!  routine ocn_init_Buttermilk_bay_destroy_topo_fields
+!
+!> \brief   Topography field cleanup routine
+!> \author  D. Wirasaet and S. Brus
+!> \date    Sep 2023
+!> \details
+!>  This routine destroys the fields that were created to hold topography
+!>  initial condition information
+!>
+!>  NOTE: adapteed from ocn_init_global_ocaen_destroy_topo_fileds
+!-----------------------------------------------------------------------
+
+    subroutine ocn_init_Buttermilk_bay_destroy_topo_fields()!{{{
+       implicit none
+
+       deallocate(topoIC % array)
+       deallocate(topoLat % array)
+       deallocate(topoLon % array)
+    end subroutine ocn_init_Buttermilk_bay_destroy_topo_fields!}}}
+
+
+!***********************************************************************
+!
+!  routine ocn_init_Buttermilk_bay
+!
+!> \brief   Validation for this initial condition
+!> \author  D. Wirasaet and Steven Brus
+!> \date   Sep 2022
+!> \details
+!>  This routine validates the configuration options for this case.
+!>
+!-----------------------------------------------------------------------
+   subroutine ocn_init_validate_Buttermilk_bay(configPool, packagePool, iocontext, iErr)!{{{
+      implicit none
+
+   !--------------------------------------------------------------------
+      type (mpas_pool_type), intent(inout) :: configPool, packagePool
+      type (mpas_io_context_type), intent(inout), target :: iocontext
+
+      integer, intent(out) :: iErr
+
+      ! character (len=StrKIND), pointer :: config_init_configuration
+      integer, pointer :: config_vert_levels, config_Buttermilk_bay_vert_levels
+      integer, pointer :: config_subgrid_table_levels, config_Buttermilk_bay_subgrid_table_levels
+
+
+      type (mpas_io_context_type), pointer :: iocontext_ptr
+      type (MPAS_IO_Handle_type) :: inputFile
+      character (len=StrKIND), pointer :: config_init_configuration, &
+                                          config_Buttermilk_bay_topography_source, &
+                                          config_Buttermilk_bay_topography_file, &
+                                          config_Buttermilk_bay_topography_nlat_dimname, &
+                                          config_Buttermilk_bay_topography_nlon_dimname
+
+      iErr = 0
+
+      call mpas_pool_get_config(configPool, 'config_init_configuration', config_init_configuration)
+
+      if(config_init_configuration .ne. trim('buttermilk_bay')) return
+
+      iocontext_ptr => iocontext
+
+      call mpas_pool_get_config(configPool, 'config_vert_levels', config_vert_levels)
+      call mpas_pool_get_config(configPool,  &
+                                'config_Buttermilk_bay_vert_levels', &
+                                config_Buttermilk_bay_vert_levels)
+
+      if(config_vert_levels <= 0 .and. config_Buttermilk_bay_vert_levels > 0) then
+         config_vert_levels = config_Buttermilk_bay_vert_levels
+      else if (config_vert_levels <= 0) then
+         call mpas_log_write( 'Validation failed for Buttermilk bay.'// &
+           'Not given a usable value for vertical levels.', MPAS_LOG_CRIT)
+         iErr = 1
+      end if
+
+      call mpas_pool_get_config(configPool, 'config_subgrid_table_levels', config_subgrid_table_levels)
+      call mpas_pool_get_config(configPool, &
+                                'config_Buttermilk_bay_subgrid_table_levels', &
+                                config_Buttermilk_bay_subgrid_table_levels)
+      if (config_subgrid_table_levels <=0 .and. config_Buttermilk_bay_subgrid_table_levels >0) then
+         config_subgrid_table_levels = config_Buttermilk_bay_subgrid_table_levels
+      else if (config_subgrid_table_levels <=0) then
+         call mpas_log_write( 'Validation failed for Buttermilk bay.'// &
+           'Not given a usable value for subgrid table levels.', MPAS_LOG_CRIT)
+         iErr = 1
+      end if
+
+
+      !---------- adapted from ocn_init_validate_global_ocean
+      call mpas_pool_get_config(configPool, 'config_Buttermilk_bay_topography_source', &
+                                config_Buttermilk_bay_topography_source)
+      call mpas_pool_get_config(configPool, 'config_Buttermilk_bay_topography_file', &
+                                config_Buttermilk_bay_topography_file)
+      call mpas_pool_get_config(configPool, 'config_Buttermilk_bay_topography_nlat_dimname', &
+                                config_Buttermilk_bay_topography_nlat_dimname)
+      call mpas_pool_get_config(configPool, 'config_Buttermilk_bay_topography_nlon_dimname', &
+                                config_Buttermilk_bay_topography_nlon_dimname)
+
+
+      call mpas_log_write( config_Buttermilk_bay_topography_source  )
+      call mpas_log_write( config_Buttermilk_bay_topography_file )
+      call mpas_log_write( config_Buttermilk_bay_topography_nlat_dimname )
+      call mpas_log_write( config_Buttermilk_bay_topography_nlon_dimname )
+
+      if (config_Buttermilk_bay_topography_source /= 'latlon_file' .and. &
+          config_Buttermilk_bay_topography_source /= 'xy_file') then
+         call mpas_log_write( 'Unexpected value for &
+                   config_Buttermilk_bay_topography_source: ' &
+             // trim(config_Buttermilk_bay_topography_source), MPAS_LOG_CRIT)
+         iErr = 1
+         return
+      end if
+
+      if (config_Buttermilk_bay_topography_file == 'none' .and. &
+          (config_Buttermilk_bay_topography_source == 'latlon_file' .or. &
+          config_Buttermilk_bay_topography_source == 'xy_file') ) then
+         call mpas_log_write( 'Validation failed for Buttermilk bay test case. ' &
+           // 'Invalid filename for config_Buttermilk_bay_topography_file ' &
+           // config_Buttermilk_bay_topography_file // '  ' &
+           // config_Buttermilk_bay_topography_source , MPAS_LOG_CRIT)
+         iErr = 1
+         return
+      end if
+
+      call mpas_log_write( ' in ocn_init_validate_buttermilk_bay '//config_Buttermilk_bay_topography_source )
+
+      if (config_Buttermilk_bay_topography_source == 'latlon_file' .or. &
+           config_Buttermilk_bay_topography_source == 'xy_file'  ) then
+
+         inputFile = MPAS_io_open(trim(config_Buttermilk_bay_topography_file), &
+                        MPAS_IO_READ, MPAS_IO_NETCDF, iocontext_ptr, ierr=iErr)
+         if (iErr /= 0) then
+            call mpas_log_write( 'could not open file '// &
+                   trim(config_Buttermilk_bay_topography_file), MPAS_LOG_CRIT)
+            return
+         end if
+
+         call MPAS_io_inq_dim(inputFile, config_Buttermilk_bay_topography_nlat_dimname, nLatTopo, iErr)
+         call MPAS_io_inq_dim(inputFile, config_Buttermilk_bay_topography_nlon_dimname, nLonTopo, iErr)
+
+         call MPAS_io_close(inputFile, iErr)
+
+      end if
+      !----------
+
+      call mpas_log_write( ' Done ocn_init_validate_buttermilk_bay ' )
+   !--------------------------------------------------------------------
+
+   end subroutine ocn_init_validate_ButterMilk_bay!}}}
+
+!***********************************************************************
+
+end module ocn_init_Buttermilk_bay
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_mode.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_mode.F
index 2ffa19a66465..0c73a0f83012 100644
--- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_mode.F
+++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_mode.F
@@ -63,8 +63,8 @@ module ocn_init_mode
    use ocn_init_mixed_layer_eddy
    use ocn_init_transport_tests
    use ocn_init_test_sht
-
    use ocn_init_parabolic_bowl
+   use ocn_init_Buttermilk_bay
 
    implicit none
    private
@@ -310,6 +310,7 @@ function ocn_init_mode_run(domain) result(iErr)!{{{
       call ocn_init_setup_transport_tests(domain, ierr)
       call ocn_init_setup_test_sht(domain, ierr)
       call ocn_init_setup_parabolic_bowl(domain, iErr)
+      call ocn_init_setup_Buttermilk_bay(domain, iErr) 
       !call ocn_init_setup_TEMPLATE(domain, ierr)
 
       call mpas_log_write( ' Completed setup of: ' // trim(config_init_configuration))
@@ -429,6 +430,8 @@ subroutine ocn_init_mode_validate_configuration(configPool, packagePool, ioconte
  
       call ocn_init_validate_parabolic_bowl(configPool, packagePool, iocontext, iErr=err_tmp)
       iErr = ior(iErr, err_tmp)
+      call ocn_init_validate_Buttermilk_bay(configPool, packagePool, iocontext, iErr=err_tmp) 
+      iErr = ior(iErr, err_tmp)
 
       ! call ocn_init_validate_TEMPLATE(configPool, packagePool, iocontext, iErr=err_tmp)
       ! iErr = ior(iErr, err_tmp)
diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F
index 43fb215ed230..a67da7e52acf 100644
--- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F
+++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_parabolic_bowl.F
@@ -33,6 +33,8 @@ module ocn_init_parabolic_bowl
    use ocn_config
    use ocn_init_vertical_grids
    use ocn_init_cell_markers
+   use ocn_subgrid
+   use ocn_init_subgrid
 
    implicit none
    private
@@ -132,6 +134,9 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{
     real (kind=RKIND):: xshift = 0.0, yshift = 0.0 
     real (kind=RKIND) :: layerThicknessEdgeAverage
     real (kind=RKIND), dimension(:,:), allocatable :: rSubgridPoints, sSubgridPoints
+    real (kind=RKIND), dimension(:), allocatable :: subgridBathymetryValues, subgridAreas
+    real (kind=RKIND), dimension(:), allocatable :: subgridSshValues
+    real (kind=RKIND), dimension(:), allocatable :: subgridUValues, subgridVValues
     real (kind=RKIND), dimension(:), allocatable :: uVelocityAverage, vVelocityAverage
     integer :: nSubgridCell, nSubgridEdge, nSubgridVertex
     integer :: nSubgridTriPerSlice
@@ -142,6 +147,12 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{
     integer :: i,j
     real (kind=RKIND) :: deltaZ
 
+    integer:: nsubgridCellEdge, iEdgeSegment
+    real (kind=RKIND), dimension(:,:), allocatable :: cellEdgeBathymetryValues
+    real (kind=RKIND), dimension(:), allocatable:: dsEdge
+    integer:: jj
+    real (kind=RKIND), dimension(:), allocatable:: xSubgridCell, ySubgridCell
+
     iErr = 0
 
     if(config_init_configuration .ne. trim('parabolic_bowl')) return
@@ -212,137 +223,180 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{
 
 
     block_ptr => domain % blocklist
-    do while(associated(block_ptr))
-       call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool)
-       call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool)
-       call mpas_pool_get_subpool(block_ptr % structs, 'verticalMesh', verticalMeshPool)
-       call mpas_pool_get_subpool(statePool, 'tracers', tracersPool)
-
-       call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) ;
-       call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) ;
-       call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve) ; 
-       call mpas_pool_get_dimension(meshPool, 'nVerticesSolve', nVerticesSolve) ; 
-       call mpas_pool_get_dimension(meshPool, 'maxEdges', maxEdges)
-      
-       call mpas_pool_get_dimension(tracersPool, 'index_temperature', index_temperature)
-       call mpas_pool_get_dimension(tracersPool, 'index_salinity', index_salinity)
+    call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool)
+    call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool)
+    call mpas_pool_get_subpool(block_ptr % structs, 'verticalMesh', verticalMeshPool)
+    call mpas_pool_get_subpool(statePool, 'tracers', tracersPool)
+
+    call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) ;
+    call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) ;
+    call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve) ; 
+    call mpas_pool_get_dimension(meshPool, 'nVerticesSolve', nVerticesSolve) ; 
+    call mpas_pool_get_dimension(meshPool, 'maxEdges', maxEdges)
+    
+    call mpas_pool_get_dimension(tracersPool, 'index_temperature', index_temperature)
+    call mpas_pool_get_dimension(tracersPool, 'index_salinity', index_salinity)
+
+    call mpas_pool_get_array(meshPool, 'xCell', xCell)
+    call mpas_pool_get_array(meshPool, 'yCell', yCell)
+    call mpas_pool_get_array(meshPool, 'refBottomDepth', refBottomDepth)
+    call mpas_pool_get_array(meshPool, 'vertCoordMovementWeights', vertCoordMovementWeights)
+    call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth)
+    call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell)
+    call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
+    call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
+    call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
+
+    call mpas_pool_get_array(meshPool, 'verticesOnCell', verticesOnCell)
+    call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge)
+    call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
+    call mpas_pool_get_array(meshPool, 'cellsOnVertex', cellsOnVertex)
+
+    call mpas_pool_get_array(meshPool, 'fCell', fCell)
+    call mpas_pool_get_array(meshPool, 'fEdge', fEdge)
+    call mpas_pool_get_array(meshPool, 'fVertex', fVertex)
+    
+    call mpas_pool_get_array(meshPool, 'xEdge', xEdge ) 
+    call mpas_pool_get_array(meshPool, 'yEdge', yEdge ) 
+    call mpas_pool_get_array(meshPool, 'xVertex', xVertex ) 
+    call mpas_pool_get_array(meshPool, 'yVertex', yVertex ) 
 
-       call mpas_pool_get_array(meshPool, 'xCell', xCell)
-       call mpas_pool_get_array(meshPool, 'yCell', yCell)
-       call mpas_pool_get_array(meshPool, 'refBottomDepth', refBottomDepth)
-       call mpas_pool_get_array(meshPool, 'vertCoordMovementWeights', vertCoordMovementWeights)
-       call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth)
-       call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell)
-       call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
-       call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
-       call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
-
-       call mpas_pool_get_array(meshPool, 'verticesOnCell', verticesOnCell)
-       call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge)
-       call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
-       call mpas_pool_get_array(meshPool, 'cellsOnVertex', cellsOnVertex)
-
-       call mpas_pool_get_array(meshPool, 'fCell', fCell)
-       call mpas_pool_get_array(meshPool, 'fEdge', fEdge)
-       call mpas_pool_get_array(meshPool, 'fVertex', fVertex)
-       
-       call mpas_pool_get_array(meshPool, 'xEdge', xEdge ) 
-       call mpas_pool_get_array(meshPool, 'yEdge', yEdge ) 
-       call mpas_pool_get_array(meshPool, 'xVertex', xVertex ) 
-       call mpas_pool_get_array(meshPool, 'yVertex', yVertex ) 
+    call mpas_pool_get_array(statePool, 'zMid', zMid, 1) ;
 
-       call mpas_pool_get_array(statePool, 'zMid', zMid, 1) ;
+    call mpas_pool_get_array(statePool, 'ssh', ssh, 1)
+    call mpas_pool_get_array(meshPool, 'edgeNormalVectors', edgeNormalVectors ) ;
+    call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity ) ;
 
-       call mpas_pool_get_array(statePool, 'ssh', ssh, 1)
-       call mpas_pool_get_array(meshPool, 'edgeNormalVectors', edgeNormalVectors ) ;
-       call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity ) ;
+    call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 1)
+    call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1)
 
-       call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 1)
-       call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1)
+    call mpas_pool_get_array(verticalMeshPool, 'refZMid', refZMid)
+    call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness)
 
-       call mpas_pool_get_array(verticalMeshPool, 'refZMid', refZMid)
-       call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness)
 
+    ! if config_parabolic_bowl_adjust_domain_center == .true.,
+    ! Adjust center of the mesh so that its center is located at (0,0)
+    if ( config_parabolic_bowl_adjust_domain_center ) then
+       xCell = xCell - xshift ; 
+       yCell = yCell - yshift ;
 
-       ! if config_parabolic_bowl_adjust_domain_center == .true.,
-       ! Adjust center of the mesh so that its center is located at (0,0)
-       if ( config_parabolic_bowl_adjust_domain_center ) then
-          xCell = xCell - xshift ; 
-          yCell = yCell - yshift ;
+       xEdge = xEdge - xshift ;
+       yEdge = yEdge - yshift ;
 
-          xEdge = xEdge - xshift ;
-          yEdge = yEdge - yshift ;
+       xVertex = xVertex - xshift ;
+       yVertex = yVertex - yshift ; 
+    end if
 
-          xVertex = xVertex - xshift ;
-          yVertex = yVertex - yshift ; 
-       end if
 
+    ! Initlialze vector
+    call mpas_initialize_vectors(meshPool) ; 
 
-       ! Initlialze edge normal vectors
-       call mpas_initialize_vectors(meshPool) ; 
 
+    minLevelCell(:) = 1
+    do iCell = 1, nCellsSolve
+      ! Set up vertical grid
+      maxLevelCell(iCell) = nVertLevels ; ! sigma coordinates 
+    end do
 
-       minLevelCell(:) = 1
-       do iCell = 1, nCellsSolve
-         ! Set up vertical grid
-         maxLevelCell(iCell) = nVertLevels ; ! sigma coordinates 
-        end do
 
+    do iCell = 1, nCellsSolve
+      
+       ! Set temperature
+       activeTracers(index_temperature, :, iCell) = 10.0_RKIND
+
+       ! Set salinity
+       activeTracers(index_salinity, :, iCell) = 30.0_RKIND
+
+      ! Set Coriolis parameters, if other than zero
+      fCell(iCell) = config_parabolic_bowl_coriolis_parameter ; 
+    end do
+
+    do iEdge = 1, nEdgesSolve
+      fEdge(iEdge) = config_parabolic_bowl_coriolis_parameter ;  
+    end do
+
+    do iVertex = 1, nVerticesSolve
+      fVertex(iVertex) =  config_parabolic_bowl_coriolis_parameter ;
+    end do
+
+    allocate(uVelocityAverage(nEdgesSolve))
+    allocate(vVelocityAverage(nEdgesSolve))
+    
+    if (config_use_subgrid_wetting_drying) then
+
+       call ocn_subgrid_init(domain,iErr)
+       call ocn_init_subgrid_calculations(domain, &
+                                          ocn_init_parabolic_bowl_bathymetry, &
+                                          ocn_init_parabolic_bowl_velocity, &
+                                          ocn_init_parabolic_bowl_ssh, &
+                                          config_parabolic_bowl_subgrid_refinement_level, &
+                                          config_parabolic_bowl_subgrid_edge_bathymetry_max_pixel, &
+                                          config_parabolic_bowl_subgrid_use_thin_layer, &
+                                          uVelocityAverage, &
+                                          vVelocityAverage, &
+                                          iErr)
+    end if
+
+    ! Find max bottom depth
+    maxBottomDepth = maxval( bottomDepth ) ; 
+    minBottomDepth = minval( bottomDepth ) ; 
+    call mpas_dmpar_max_real( domain % dminfo, maxBottomDepth, globalMaxBottomDepth ) ;
+    call mpas_dmpar_min_real( domain % dminfo, minBottomDepth, globalMinBottomDepth ) ;
+
+    ! Set refBottomDepth and refZMid
+    do k = 1, nVertLevels
+      refBottomDepth(k) = globalMaxBottomDepth*interfaceLocations(k+1) ; 
+      refZMid(k) = -0.5_RKIND*( interfaceLocations(k+1) + interfaceLocations(k))*globalMaxBottomDepth ; 
+    end do
+
+    ! Set vertCoordMovementWeights
+    vertCoordMovementWeights(:) = 1.0_RKIND
+
+    ! Set velocity
+    do iEdge = 1, nEdgesSolve
+
+      if (config_use_subgrid_wetting_drying) then
+         do k = 1, nVertLevels 
+            normalVelocity(k,iEdge) = uVelocityAverage(iEdge)*edgeNormalVectors(1,iEdge) &
+                                    + vVelocityAverage(iEdge)*edgeNormalVectors(2,iEdge) ;
+         end do
+      else
+         call ocn_init_parabolic_bowl_velocity(xEdge(iEdge), yEdge(iEdge), uu, vv)
+         do k = 1, nVertLevels 
+            normalVelocity(k,iEdge) = uu*edgeNormalVectors(1,iEdge) + vv*edgeNormalVectors(2,iEdge) ;
+         end do
+      end if
+    end do
+
+    ! Set layer thickness and ssh
+    if (config_use_wetting_drying) then
+
+      do iCell = 1, nCellsSolve
+        ! Set up vertical grid
+        maxLevelCell(iCell) = nVertLevels ; ! sigma coordinates 
+       end do
 
        do iCell = 1, nCellsSolve
-         
-          ! Set temperature
-          activeTracers(index_temperature, :, iCell) = 10.0_RKIND
-
-          ! Set salinity
-          activeTracers(index_salinity, :, iCell) = 30.0_RKIND
-
-         ! Set Coriolis parameters
-         fCell(iCell) = config_parabolic_bowl_coriolis_parameter ; 
-      end do
-
-      do iEdge = 1, nEdgesSolve
-         fEdge(iEdge) = config_parabolic_bowl_coriolis_parameter ;  
-      end do
-
-      do iVertex = 1, nVerticesSolve
-         fVertex(iVertex) =  config_parabolic_bowl_coriolis_parameter ;
-      end do
-
-
-      ! Find max bottom depth
-      maxBottomDepth = maxval( bottomDepth ) ; 
-      minBottomDepth = minval( bottomDepth ) ; 
-      call mpas_dmpar_max_real( domain % dminfo, maxBottomDepth, globalMaxBottomDepth ) ;
-      call mpas_dmpar_min_real( domain % dminfo, minBottomDepth, globalMinBottomDepth ) ;
-
-      ! Set refBottomDepth and refZMid
-      do k = 1, nVertLevels
-        refBottomDepth(k) = globalMaxBottomDepth*interfaceLocations(k+1) ; 
-        refZMid(k) = -0.5_RKIND*( interfaceLocations(k+1) + interfaceLocations(k))*globalMaxBottomDepth ; 
-      end do
-
-      ! Set vertCoordMovementWeights
-      vertCoordMovementWeights(:) = 1.0_RKIND
-
-      ! Set velocity
-      do iEdge = 1, nEdgesSolve
-
-        call ocn_init_parabolic_bowl_velocity(xEdge(iEdge), yEdge(iEdge), uu, vv)
-        do k = 1, nVertLevels 
-           normalVelocity(k,iEdge) = uu*edgeNormalVectors(1,iEdge) + vv*edgeNormalVectors(2,iEdge) ;
-        end do
-      end do
-
-      ! Set layer thickness and ssh
-      if (config_use_wetting_drying) then
-
-         do iCell = 1, nCellsSolve
-           !
-           ! make sure depth is thick enough via ssh = TOTAL_DEPTH - bottomDepth
-           ! add a thin layer of nlayer*config_drying_min_cellhight 
-           ! 
+         !
+         ! make sure depth is thick enough via ssh = TOTAL_DEPTH - bottomDepth
+         ! add a thin layer of nlayer*config_drying_min_cellhight 
+         ! 
  
+         if (config_use_subgrid_wetting_drying) then
+
+           call ocn_subgrid_ssh_lookup(layerThickness(1,iCell),&
+                                       subgridWetVolumeCellTable(:,iCell),&
+                                       subgridSshCellTableRange(:,iCell),&
+                                       bottomDepth(iCell),&
+                                       subgridCellBathymetryMin(iCell),&
+                                       ssh(iCell))
+           !call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & 
+           !                                subgridWetVolumeCellTable(:,iCell), &
+           !                                subgridSshCellTableRange(:,iCell),&
+           !                                bottomDepth(iCell),&
+           !                                LayerThickness(1,iCell))
+
+         else
            call ocn_init_parabolic_bowl_bathymetry(xCell(iCell),yCell(iCell),bottomDepth(iCell))
            call ocn_init_parabolic_bowl_ssh(xCell(iCell),yCell(iCell),bottomDepth(iCell),ssh(iCell))
            ssh(iCell) = - bottomDepth(iCell) + &
@@ -359,17 +413,25 @@ subroutine ocn_init_setup_parabolic_bowl(domain, iErr)!{{{
                  realArgs=(/layerThickness(k,iCell)/))
              end if
            end do
-         
-         
-           do k = 1, maxLevelCell(iCell)
-             restingThickness(k,iCell) = bottomDepth(iCell)/maxLevelCell(iCell)
-           end do
+         endif
+       
+       
+         do k = 1, maxLevelCell(iCell)
+           restingThickness(k,iCell) = bottomDepth(iCell)/maxLevelCell(iCell)
          end do
-         
-      end if 
-
-      block_ptr => block_ptr % next
-    end do
+       end do
+       
+    end if 
+
+    if (config_use_subgrid_wetting_drying) then
+    do iCell = 1,nCellsSolve
+      call ocn_subgrid_layer_thickness_lookup(ssh(iCell), & 
+                                              subgridWetVolumeCellTable(:,iCell), &
+                                              subgridSshCellTableRange(:,iCell),&
+                                              bottomDepth(iCell),&
+                                              subgridLayerThicknessDebug(iCell))
+    enddo
+    endif
 
     deallocate(interfaceLocations)
     !--------------------------------------------------------------------
@@ -436,6 +498,7 @@ subroutine ocn_init_parabolic_bowl_ssh(x, y, bottomDepth, ssh)!{{{
  
       ssh = (sqrtOneMC2/oneMC) - 1.0_RKIND - ((RR**2)/(LL**2))*( (oneMC2/(oneMC**2)) - 1.0_RKIND ) ;
       ssh = config_parabolic_bowl_b0*ssh ; 
+      !ssh = - bottomDepth + max(ssh + bottomDepth, config_drying_min_cell_height + eps)
 
    !--------------------------------------------------------------------
 
@@ -506,6 +569,7 @@ subroutine ocn_init_validate_parabolic_bowl(configPool, packagePool, iocontext,
 
       character (len=StrKIND), pointer :: config_init_configuration
       integer, pointer :: config_vert_levels, config_parabolic_bowl_vert_levels
+      integer, pointer :: config_subgrid_table_levels, config_parabolic_bowl_subgrid_table_levels
    
       iErr = 0
 
@@ -519,7 +583,17 @@ subroutine ocn_init_validate_parabolic_bowl(configPool, packagePool, iocontext,
       if(config_vert_levels <= 0 .and. config_parabolic_bowl_vert_levels > 0) then
          config_vert_levels = config_parabolic_bowl_vert_levels
       else if (config_vert_levels <= 0) then
-         call mpas_log_write( 'Validation failed for para_bowl. Not given a usable value for vertical levels.', MPAS_LOG_CRIT)
+         call mpas_log_write( 'Validation failed for parabolic_bowl. Not given a usable value for vertical levels.', MPAS_LOG_CRIT)
+         iErr = 1
+      end if
+
+      call mpas_pool_get_config(configPool, 'config_subgrid_table_levels', config_subgrid_table_levels)
+      call mpas_pool_get_config(configPool, 'config_parabolic_bowl_subgrid_table_levels', config_parabolic_bowl_subgrid_table_levels)
+
+      if (config_subgrid_table_levels <= 0 .and. config_parabolic_bowl_subgrid_table_levels > 0) then
+         config_subgrid_table_levels = config_parabolic_bowl_subgrid_table_levels  
+      else if (config_subgrid_table_levels <= 0) then
+         call mpas_log_write( 'Validation failed for parabolic_bowl. Not given a usable value for subgrid table levels.', MPAS_LOG_CRIT)
          iErr = 1
       end if
 
diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_subgrid.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_subgrid.F
new file mode 100644
index 000000000000..37c3e2a273d4
--- /dev/null
+++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_subgrid.F
@@ -0,0 +1,941 @@
+! Copyright (c) 2013,  Los Alamos National Security, LLC (LANS)
+! and the University Corporation for Atmospheric Research (UCAR).
+!
+! Unless noted otherwise source code is licensed under the BSD license.
+! Additional copyright and license information can be found in the LICENSE file
+! distributed with this code, or at http://mpas-dev.github.com/license.html
+!
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  ocn_init_subgrid
+!
+!> \brief  Calculate the subgrid information 
+!> \author D. Wirasaet, S. Brus
+!> \date   May-June 2022
+!> \details
+!>  This module contains the routines for calculating the information
+!>  needed to use subgrid corrctions within a forward model run.
+!>
+!-----------------------------------------------------------------------
+
+module ocn_init_subgrid
+
+   use mpas_kind_types
+   use mpas_io_units
+   use mpas_derived_types
+   use mpas_pool_routines
+   use mpas_constants
+   use mpas_dmpar
+
+   use ocn_constants
+   use ocn_config
+   use ocn_subgrid
+
+   use mpas_constants
+   use mpas_io
+   use mpas_io_streams
+   use mpas_stream_manager
+
+
+   implicit none
+   private
+   save
+
+   !--------------------------------------------------------------------
+   !
+   ! Public parameters
+   !
+   !--------------------------------------------------------------------
+
+   !--------------------------------------------------------------------
+   !
+   ! Public member functions
+   !
+   !--------------------------------------------------------------------
+
+   public :: ocn_init_subgrid_calculations 
+
+   !--------------------------------------------------------------------
+   !
+   ! Private module variables
+   !
+   !--------------------------------------------------------------------
+
+   real(kind=RKIND), parameter:: eps = 1.0e-10_RKIND ;   
+   real(kind=RKIND) :: subgrid_thin_layer
+
+   abstract interface
+      subroutine bathymetry_function(x, y, b)
+         use mpas_kind_types, only: rkind
+         implicit none
+         real (kind=RKIND), intent(in) :: x
+         real (kind=RKIND), intent(in) :: y
+         real (kind=RKIND), intent(out) :: b
+      end subroutine
+      subroutine velocity_function(x, y, u, v)
+         use mpas_kind_types, only: rkind
+         implicit none
+         real (kind=RKIND), intent(in) :: x
+         real (kind=RKIND), intent(in) :: y
+         real (kind=RKIND), intent(out) :: u
+         real (kind=RKIND), intent(out) :: v
+      end subroutine
+      subroutine ssh_function(x, y, b, z)
+         use mpas_kind_types, only: rkind
+         implicit none
+         real (kind=RKIND), intent(in) :: x
+         real (kind=RKIND), intent(in) :: y
+         real (kind=RKIND), intent(in) :: b
+         real (kind=RKIND), intent(out) :: z
+      end subroutine
+   end interface
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+!  routine ocn_init_subgrid_calculations
+!
+!> \brief   Calculations for subgrid look-up tables 
+!> \author  D. Wirasaet and S. Brus
+!> \date    May-June 2022
+!> \details
+!>    This routine performs the vertical integration of the wet fraction
+!>    over the cell, edge, and vertex control volumes to create
+!>    lookup tables for the relationship between ssh and wet volume per
+!>    unit area.
+!
+!-----------------------------------------------------------------------
+
+  subroutine ocn_init_subgrid_calculations(domain, problem_bathymetry, problem_velocity, problem_ssh, &
+                                           subgrid_refinement_level, &
+                                           subgrid_edge_bathymetry_max_pixel, &
+                                           subgrid_use_thin_layer, &
+                                           uVelocityAverage, vVelocityAverage, iErr)!{{{
+    
+    implicit none
+    !--------------------------------------------------------------------
+
+    type (domain_type), intent(inout) :: domain
+    integer, intent(out) :: iErr
+    procedure(bathymetry_function) :: problem_bathymetry
+    procedure(velocity_function) :: problem_velocity
+    procedure(ssh_function) :: problem_ssh
+    integer, intent(in) :: subgrid_refinement_level
+    logical, intent(in) :: subgrid_edge_bathymetry_max_pixel
+    logical, intent(in) :: subgrid_use_thin_layer
+    
+    real (kind=RKIND), dimension(:), intent(out) :: uVelocityAverage, vVelocityAverage
+
+    type (block_type), pointer :: block_ptr
+    type (mpas_pool_type), pointer :: meshPool
+    type (mpas_pool_type), pointer :: statePool
+
+    ! local variables
+    integer :: iCell, iEdge, iVertex, k, idx
+
+    ! Define dimension pointers
+    integer, pointer :: nCellsSolve, nEdgesSolve, nVerticesSolve 
+    integer, pointer :: maxEdges
+
+    ! Define variable pointers
+    integer, dimension(:), pointer :: nEdgesOnCell
+    integer, dimension(:,:), pointer :: verticesOnCell, verticesOnEdge
+    integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex
+    real (kind=RKIND), dimension(:), pointer :: xCell, yCell, bottomDepth 
+    real (kind=RKIND), dimension(:,:), pointer:: edgeNormalVectors
+    real (kind=RKIND), dimension(:), pointer:: xEdge, yEdge, xVertex, yVertex
+    
+    real (kind=RKIND), dimension(:), pointer :: ssh
+    real (kind=RKIND), dimension(:,:), pointer :: normalVelocity
+    real (kind=RKIND), dimension(:,:), pointer :: layerThickness
+    real (kind=RKIND), dimension(:), pointer :: areaCell
+
+    real (kind=RKIND) :: layerThicknessEdgeAverage
+    real (kind=RKIND), dimension(:,:), allocatable :: rSubgridPoints, sSubgridPoints
+    real (kind=RKIND), dimension(:), allocatable :: subgridBathymetryValues, subgridAreas
+    real (kind=RKIND), dimension(:), allocatable :: subgridSshValues
+    real (kind=RKIND), dimension(:), allocatable :: subgridUValues, subgridVValues
+    integer :: nSubgridCell, nSubgridEdge, nSubgridVertex
+    integer :: nSubgridTriPerSlice
+    integer :: v1, v2
+    integer :: c1, c2
+    real (kind=RKIND) :: x(3), y(3)
+    integer :: slice, nSlice
+    real (kind=RKIND) :: deltaZ
+
+
+
+    integer:: i, j, jj
+    integer:: nsubgridCellEdge, iEdgeSegment
+    real (kind=RKIND):: pi
+    real (kind=RKIND), dimension(:,:), allocatable :: cellEdgeBathymetryValues
+    real (kind=RKIND), dimension(:), allocatable:: dsEdge
+    real (kind=RKIND), dimension(:), allocatable:: xSubgridCell, ySubgridCell
+    real (kind=RKIND):: bathymetryMin, bathymetryMax 
+    iErr = 0
+
+
+    block_ptr => domain % blocklist
+    call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool)
+    call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool)
+
+    call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) ;
+    call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve) ; 
+    call mpas_pool_get_dimension(meshPool, 'nVerticesSolve', nVerticesSolve) ; 
+    call mpas_pool_get_dimension(meshPool, 'maxEdges', maxEdges)
+    
+    call mpas_pool_get_array(meshPool, 'xCell', xCell)
+    call mpas_pool_get_array(meshPool, 'yCell', yCell)
+    call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth)
+    call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
+    call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
+
+    call mpas_pool_get_array(meshPool, 'verticesOnCell', verticesOnCell)
+    call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge)
+    call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
+    call mpas_pool_get_array(meshPool, 'cellsOnVertex', cellsOnVertex)
+
+    call mpas_pool_get_array(meshPool, 'xEdge', xEdge ) 
+    call mpas_pool_get_array(meshPool, 'yEdge', yEdge ) 
+    call mpas_pool_get_array(meshPool, 'xVertex', xVertex ) 
+    call mpas_pool_get_array(meshPool, 'yVertex', yVertex ) 
+
+    call mpas_pool_get_array(statePool, 'ssh', ssh, 1)
+    call mpas_pool_get_array(meshPool, 'edgeNormalVectors', edgeNormalVectors ) ; 
+    call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity ) ; 
+    call mpas_pool_get_array(statePool, 'layerThickness', layerThickness) ; 
+
+    if (subgrid_use_thin_layer) then
+       subgrid_thin_layer = config_drying_min_cell_height + eps
+    else
+       subgrid_thin_layer = 0.0_RKIND
+    end if
+
+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    ! Cells
+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    print*, "Begin cells", nCellsSolve
+
+    nSubgridTriPerSlice = subgrid_refinement_level**2
+    allocate(rSubgridPoints(3,maxEdges*nSubgridTriPerSlice), sSubgridPoints(3,maxEdges*nSubgridTriPerSlice))
+    call ocn_init_define_subgrid_points(subgrid_refinement_level, rSubgridPoints, sSubgridPoints)
+    allocate(subgridBathymetryValues(maxEdges*nSubgridTriPerSlice), subgridAreas(maxEdges*nSubgridTriPerSlice))
+    allocate(subgridSshValues(maxEdges*nSubgridTriPerSlice))
+    allocate(subgridUValues(maxEdges*nSubgridTriPerSlice))
+    allocate(subgridVValues(maxEdges*nSubgridTriPerSlice))
+
+    allocate(xSubgridCell(maxEdges*nSubgridTriPerSlice)) 
+    allocate(ySubgridCell(maxEdges*nSubgridTriPErSlice)) 
+    xSubgridCell = 0.0_RKIND  
+    ySubgridCell = 0.0_RKIND 
+
+    do iCell = 1,nCellsSolve
+      !
+      ! Evaluate subgrid bathymetry at centers of sub-triangles for cell slices 
+      ! (all subdivided triangles for each cell slice are gathered into 
+      !  subgridBathymetryValues and subgridAreas)
+      !---------------------------------------------------------------
+
+      nSubgridCell = 0  ! Counter for all subgrid triangles over cell slices
+      do slice = 1,nEdgesOnCell(iCell) ! Loop over cell slices
+
+        v1 = verticesOnCell(slice,iCell)
+        if (slice+1 <= nEdgesOnCell(iCell)) then
+          v2 = verticesOnCell(slice+1,iCell)
+        else
+          v2 = verticesOnCell(1,iCell)
+        endif
+
+        ! Cell slice coordinates
+        x(1) = xCell(iCell)
+        y(1) = yCell(iCell)
+
+        x(2) = xVertex(v1)
+        y(2) = yVertex(v1)
+
+        x(3) = xVertex(v2)
+        y(3) = yVertex(v2)
+
+        call ocn_init_evaluate_subgrid_data(x, y, nSubgridTriPerSlice, nSubgridCell, rSubgridPoints, sSubgridPoints, &
+                                            problem_bathymetry, problem_velocity, problem_ssh, &  
+                                            subgridBathymetryValues, subgridAreas, subgridSshValues, &
+                                            subgridUValues, subgridVValues, xSubgridCell, ysubgridCell) !{{{
+      enddo
+    
+      ! Evaluate bounds of look-up table range
+      !---------------------------------------------------------------        
+      bathymetryMin = maxval( subgridBathymetryValues(1:nSubgridCell) ) 
+      bathymetryMax = minval( subgridBathymetryValues(1:nSubgridCell) )
+
+      if ( abs(bathymetryMin - bathymetryMax) > 100.0*eps ) then
+          subgridSshCellTableRange(1,iCell) = -maxval(subgridBathymetryValues(1:nSubgridCell)) + subgrid_thin_layer
+          subgridSshCellTableRange(2,iCell) = -minval(subgridBathymetryValues(1:nSubgridCell))
+      else
+          ! flat bathy !
+          subgridSshCellTableRange(1,iCell) =  -bathymetryMin + subgrid_thin_layer
+          subgridSshCellTableRange(2,iCell) =  -bathymetryMin + 2.0*config_drying_min_cell_height  
+      endif 
+
+      ! Evaluate subgrid bathymetry
+      !---------------------------------------------------------------
+      bottomDepth(iCell) = sum(subgridBathymetryValues(1:nSubgridCell)*subgridAreas(1:nSubgridCell))/sum(subgridAreas(1:nSubgridCell))
+      subgridCellBathymetryMin(iCell) = maxval(subgridBathymetryValues(1:nSubgridCell))
+
+      ! Vertical integration of wet fraction
+      !---------------------------------------------------------------
+
+      call ocn_init_vertical_integration(iCell,subgridSshCellTableRange, nSubgridCell, subgridBathymetryValues, subgridAreas, &
+                                         subgridWetVolumeCellTable, subgridWetFractionCellTable)
+
+      ! Evaluate wet layerThickness average
+      !---------------------------------------------------------------
+        call ocn_init_grid_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThickness(1,iCell))
+!        call ocn_init_wet_average_ssh(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, ssh(iCell))!{{{
+    enddo
+
+    do iCell = 1, nCellsSolve
+        call ocn_subgrid_ssh_lookup( config_drying_min_cell_height, & 
+                                     subgridWetVolumeCellTable(:,iCell), &
+                                     subgridSshCellTableRange(:,iCell),&
+                                     bottomDepth(iCell),&
+                                     subgridCellBathymetryMin(iCell),&
+                                     subgridSShCellTableRange(3,iCell) )
+
+       if ( ssh(iCell) < subgridSshCellTableRange(3,iCell) ) then
+          ssh(iCell) = subgridSshCellTableRange(3,iCell) ;  
+       end if 
+    end do
+
+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    ! Edges
+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    print*, "Begin edges", nEdgesSolve
+
+    nSubgridCellEdge=subgrid_refinement_level 
+    allocate( cellEdgeBathymetryValues(3,nSubgridCellEdge) )   
+    allocate( dsEdge(nSubgridCellEdge) ) 
+    cellEdgeBathymetryValues = -99999 ; 
+
+
+    do iEdge = 1,nEdgesSolve
+      !--------------------------------------------------------------
+      ! Evaluate subgrid bathymetry at centers of sub-triangles for edge slices
+      ! (all subdivided triangles for each edge slice are gathered into 
+      !  subgridBathymetryValues and subgridAreas)
+      !---------------------------------------------------------------
+
+      nSlice = 0
+      do slice = 1,2
+        if (cellsOnEdge(slice,iEdge) <= nCellsSolve) then
+          nSlice = nSlice + 1
+        endif
+      enddo
+
+      nSubgridEdge = 0  ! Counter for all subgrid triangles over edge slices
+      do slice = 1,nSlice ! Loop over edge slices
+
+        ! Edge slice coordinates
+        x(1) = xVertex(verticesOnEdge(1,iEdge))
+        y(1) = yVertex(verticesOnEdge(1,iEdge))
+
+        x(2) = xVertex(verticesOnEdge(2,iEdge))
+        y(2) = yVertex(verticesOnEdge(2,iEdge))
+
+        x(3) = xCell(cellsOnEdge(slice,iEdge))
+        y(3) = yCell(cellsOnEdge(slice,iEdge))
+
+
+        call ocn_init_evaluate_subgrid_data(x, y, nSubgridTriPerSlice, nSubgridEdge, rSubgridPoints, sSubgridPoints, &
+                                            problem_bathymetry, problem_velocity, problem_ssh, &  
+                                            subgridBathymetryValues, subgridAreas, subgridSshValues, subgridUValues, subgridVValues, &
+                                            xSubgridCell, ySubgridCell)
+  
+      enddo
+
+      ! Evaluate velocity average
+      !---------------------------------------------------------------
+      call ocn_init_grid_average(nSubgridCell, subgridBathymetryValues, subgridSshValues, subgridAreas, layerThicknessEdgeAverage, &
+                                subgridUValues, subgridVValues, uVelocityAverage(iEdge), vVelocityAverage(iEdge))
+
+      if ( .NOT. subgrid_edge_bathymetry_max_pixel) then
+         ! Evaluate bounds of look-up table range
+         !---------------------------------------------------------------
+         subgridSshEdgeTableRange(1,iEdge) = -maxval(subgridBathymetryValues(1:nSubgridEdge)) + subgrid_thin_layer
+         subgridSshEdgeTableRange(2,iEdge) = -minval(subgridBathymetryValues(1:nSubgridEdge))
+
+         ! Evaluate subgrid bathymetry
+         !---------------------------------------------------------------
+         subgridEdgeBathymetryMean(iEdge) = sum(subgridBathymetryValues(1:nSubgridEdge))/real(nSubgridEdge,RKIND) 
+         subgridEdgeBathymetryMin(iEdge)  = maxval(subgridBathymetryValues(1:nSubgridEdge))
+
+         ! Vertical integration of wet fraction
+         !---------------------------------------------------------------
+         call ocn_init_vertical_integration(iEdge,subgridSshEdgeTableRange, nSubgridEdge, subgridBathymetryValues, subgridAreas, &
+                                         subgridWetVolumeEdgeTable, subgridWetFractionEdgeTable)
+      else
+         ! DW: Use the higher values of the pair of subcells along the
+         !     cell edge
+         cellEdgeBathymetryValues(1,:) = subgridBathymetryValues(1:2*nSubgridCellEdge - 1:2) 
+         if ( nslice > 1 ) then
+           cellEdgeBathymetryValues(2,:) = &
+              subgridBathymetryValues(nSubgridTriPerSlice+2*nsubgridCellEdge - 1:nSubgridTriPerSlice+1:-2)
+         else
+           cellEdgeBathymetryValues(2,:) = cellEdgeBathymetryValues(1,:)     
+         endif
+
+         do iEdgeSegment = 1, nSubgridCellEdge
+            cellEdgeBathymetryValues(3,iEdgeSegment) = minval(cellEdgeBathymetryValues(1:2,iEdgeSegment) ) 
+         end do 
+         dsEdge(:) = sqrt( (x(2) - x(1))*(x(2) - x(1)) + (y(2) - y(1))*(y(2)- y(1)) )/nSubgridCellEdge 
+
+
+         ! Evaluate bounds of look-up table range
+         !---------------------------------------------------------------
+         subgridSshEdgeTableRange(1,iEdge) = -maxval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge)) + subgrid_thin_layer
+         subgridSshEdgeTableRange(2,iEdge) = -minval(CellEdgeBathymetryValues(3,1:nSubgridCellEdge))
+
+         ! Evaluate bounds of look-up table range
+         !---------------------------------------------------------------
+         subgridEdgeBathymetryMean(iEdge) = sum(cellEdgeBathymetryValues(3,1:nSubgridCellEdge))/real(nSubgridCellEdge,RKIND) 
+         subgridEdgeBathymetryMin(iEdge)  = maxval(cellEdgeBathymetryValues(3,1:nSubgridCellEdge))
+
+         ! Vertical integration of wet fraction
+         !---------------------------------------------------------------
+         call ocn_init_vertical_integration( iEdge, subgridSshEdgeTableRange, &
+                    nSubgridCellEdge, cellEdgeBathymetryValues(3,:), dsEdge, subgridWetVolumeEdgeTable, subgridWetFractionEdgeTable )
+      endif
+
+    end do
+
+    ! find an ssh value corresponding to drying_min_cell-height 
+    ! of each edge 
+    do iEdge = 1,nEdgesSolve
+        call ocn_subgrid_ssh_lookup( config_drying_min_cell_height, & 
+                                     subgridWetVolumeEdgeTable(:,iEdge),&
+                                     subgridSshEdgeTableRange(:,iEdge),&
+                                     subgridEdgeBathymetryMean(iEdge),&
+                                     subgridEdgeBathymetryMin(iEdge),&
+                                     subgridSshEdgeTableRange(3,iEdge) )
+    end do
+
+
+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    ! Vertex
+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    print*, "Begin vertex", nVerticesSolve
+
+vertex: do iVertex = 1,nVerticesSolve
+
+      ! Evaluate subgrid bathymetry at centers of sub-triangles for vertex triangle
+      !---------------------------------------------------------------
+
+      nSlice = 0
+      do slice = 1,3
+        if (cellsOnVertex(slice,iVertex) <= nCellsSolve) then
+          nSlice = nSlice + 1
+        endif
+      enddo
+
+      if (nSlice < 3) then
+         cycle vertex
+      endif
+
+      nSubgridVertex = 0  ! Counter for all subgrid triangles over edge slices
+      do slice = 1,nSlice
+
+        c1 = cellsOnVertex(slice,iVertex)
+        if (slice < 3) then
+          c2 = cellsOnVertex(slice+1,iVertex)  
+        else
+          c2 = cellsOnVertex(1,iVertex)
+        endif
+
+        ! Vertex slice coordinates
+        x(1) = xCell(c1)
+        y(1) = yCell(c1)
+
+        x(2) = xCell(c2)
+        y(2) = yCell(c2)
+
+        x(3) = xVertex(iVertex)
+        y(3) = yVertex(iVertex)
+
+        call ocn_init_evaluate_subgrid_data(x, y, nSubgridTriPerSlice, nSubgridVertex, rSubgridPoints, sSubgridPoints, &
+                                            problem_bathymetry, problem_velocity, problem_ssh, &  
+                                            subgridBathymetryValues, subgridAreas)
+      enddo
+
+      ! Evaluate bounds of look-up table range
+      !---------------------------------------------------------------
+      subgridSshVertexTableRange(1,iVertex) = -maxval(subgridBathymetryValues(1:nSubgridVertex)) + subgrid_thin_layer
+      subgridSshVertexTableRange(2,iVertex) = -minval(subgridBathymetryValues(1:nSubgridVertex))
+
+      ! Evaluate subgrid bathymetry
+      !---------------------------------------------------------------
+      subgridVertexBathymetryMean(iVertex) = sum(subgridBathymetryValues(1:nSubgridVertex)*subgridAreas(1:nSubgridVertex))/sum(subgridAreas(1:nSubgridVertex))
+      subgridVertexBathymetryMin(iVertex)  = maxval(subgridBathymetryValues(1:nSubgridVertex))
+
+      ! Vertical integration of wet fraction
+      !---------------------------------------------------------------
+      call ocn_init_vertical_integration(iVertex,subgridSshVertexTableRange, nSubgridVertex, subgridBathymetryValues, subgridAreas, &
+                                         subgridWetVolumeVertexTable, subgridWetFractionVertexTable)
+
+    enddo vertex
+
+    ! find an ssh value corresponding to drying_min_cell-height 
+    ! of each edge 
+    do iVertex = 1,nVerticesSolve
+
+        call ocn_subgrid_ssh_lookup( config_drying_min_cell_height, & 
+                                     subgridWetVolumeVertexTable(:,iVertex),&
+                                     subgridSshVertexTableRange(:,iVertex),&
+                                     subgridVertexBathymetryMean(iVertex),&
+                                     subgridVertexBathymetryMin(iVertex),&
+                                     subgridSshVertexTableRange(3,iVertex) )
+    end do
+
+    return ;  
+  end subroutine ocn_init_subgrid_calculations!}}}
+
+!***********************************************************************
+!
+!  routine ocn_init_define_subgrid_points
+!
+!> \brief   Define subgrid points on reference triangle
+!> \author  Steven Brus 
+!> \date    November 2022 
+!> \details Gives the r and s coordinates of the subgrid triangle vertices
+!>  
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_init_define_subgrid_points(nSubgridLevel, rSubgridPoints, sSubgridPoints)!{{{
+           
+      implicit none 
+
+      integer, intent(in) :: nSubgridLevel
+      real (kind=RKIND), dimension(:,:), intent(inout) :: rSubgridPoints, sSubgridPoints
+
+      integer :: i,j,k
+      integer :: nSubgridTri
+      real (kind=RKIND) :: dx
+      real (kind=RKIND), dimension(:), allocatable :: xPoints
+
+   !--------------------------------------------------------------------
+
+      ! Equi-spaced nodes on -1,+1
+      allocate(xPoints(nSubgridLevel+1))
+      dx = 2.0_RKIND/real(nSubgridLevel,RKIND)
+      xPoints(1) = -1.0_RKIND
+      do i = 2,nSubgridLevel+1
+        xPoints(i) = xPoints(i-1)+dx   
+      enddo      
+
+      !   * 
+      !   |\
+      !   | \
+      !   |  \
+      !   |   \
+      ! s *----*
+      !   |\ u |\
+      !   | \  | \
+      !   |  \ |  \
+      !   | l \| l \
+      !   *----*----*
+      !        r
+
+      ! Trianglulate tensor product of equi-spaced nodes on reference triangle
+      ! r - horizontal coordinate (i index), s - vertical coordinate (j index)
+      nSubgridTri = 0
+      do j = 1,nSubgridLevel
+        do i = 1,nSubgridLevel+1 - j 
+
+          nSubgridTri = nSubgridTri+1 
+
+          ! lower triangle in pair (see l in triangle above)
+          rSubgridPoints(1,nSubgridTri) = xPoints(i)
+          sSubgridPoints(1,nSubgridTri) = xPoints(j)
+          
+          rSubgridPoints(2,nSubgridTri) = xPoints(i+1)
+          sSubgridPoints(2,nSubgridTri) = xPoints(j)
+
+          rSubgridPoints(3,nSubgridTri) = xPoints(i)
+          sSubgridPoints(3,nSubgridTri) = xPoints(j+1)
+
+          ! upper triangle in pair (see u in triangle above, doesn't occur next to hypotenuse)          
+          if (i < nSubgridLevel+1 - j) then
+            
+            nSubgridTri = nSubgridTri + 1
+
+            rSubgridPoints(1,nSubgridTri) = xPoints(i+1)
+            sSubgridPoints(1,nSubgridTri) = xPoints(j)
+            
+            rSubgridPoints(2,nSubgridTri) = xPoints(i+1)
+            sSubgridPoints(2,nSubgridTri) = xPoints(j+1)
+
+            rSubgridPoints(3,nSubgridTri) = xPoints(i)
+            sSubgridPoints(3,nSubgridTri) = xPoints(j+1)
+
+          endif
+
+        enddo
+      enddo
+
+      print*, "exit ocn_init_define_subgrid_points"
+
+   
+   !--------------------------------------------------------------------
+
+   end subroutine ocn_init_define_subgrid_points!}}}
+
+!***********************************************************************
+!
+!  routine ocn_init_evaluate_subgrid_data
+!
+!> \brief   Evaluate subgrid infromation 
+!> \author  Steven Brus
+!> \date    November 2022
+!> \details Evaluate data on subgrid triangles for a given triangluar grid cell region 
+!>  
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_init_evaluate_subgrid_data(xTri, yTri, nSubgridTri, nSubgridCV, rSubgridPoints, sSubgridPoints, & 
+                                             problem_bathymetry, problem_velocity, problem_ssh, &
+                                             subgridBathymetryValues, subgridAreas, &
+                                             subgridSshValues, subgridUValues, subgridVValues, &
+                                             xSubgridCell, ySubgridCell) !{{{
+
+           
+      implicit none 
+      real (kind=RKIND), intent(inout) :: xTri(3), yTri(3)
+      integer, intent(in) :: nSubgridTri
+      integer, intent(inout) :: nSubgridCV
+
+      procedure(bathymetry_function) :: problem_bathymetry
+      procedure(velocity_function) :: problem_velocity
+      procedure(ssh_function) :: problem_ssh
+
+      real (kind=RKIND), dimension(:,:), intent(in) :: rSubgridPoints, sSubgridPoints
+      real (kind=RKIND), dimension(:), intent(inout) :: subgridBathymetryValues
+      real (kind=RKIND), dimension(:), intent(inout) :: subgridAreas
+      real (kind=RKIND), dimension(:), intent(inout),optional :: subgridSshValues
+      real (kind=RKIND), dimension(:), intent(inout),optional :: subgridUValues, subgridVValues
+
+      real (kind=RKIND), dimension(:), intent(inout), optional :: xSubgridCell, ySubgridCell
+
+      real (kind=RKIND) :: rCenter, sCenter
+      real (kind=RKIND) :: xCenter, yCenter
+      real (kind=RKIND) :: xSubgridPoints(3), ySubgridPoints(3)
+      real (kind=RKIND) :: x, y
+      real (kind=RKIND) :: area
+
+
+      integer :: iPt,i
+
+   !--------------------------------------------------------------------
+
+     ! Coordinates of physical triangle
+     ! (ensure counter-clockwise numering)
+     call ocn_init_tri_area(xTri(:), yTri(:), area)
+     if (area < 0.0_RKIND) then
+      x = xTri(1)
+      y = yTri(1) 
+      xTri(1) = xTri(2) 
+      yTri(1) = yTri(2)
+      xTri(2) = x ;
+      yTri(2) = y ; 
+      area = abs(area) ; 
+     endif
+
+      do iPt = 1,nSubgridTri
+
+        ! Counter over all subcells within cell/edge/vertex control volume
+        nSubgridCV = nSubgridCV + 1
+
+        ! Center sub-triangle (on reference triangle)
+        rCenter = sum(rSubgridPoints(:,iPt))/3.0_RKIND
+        sCenter = sum(sSubgridPoints(:,iPt))/3.0_RKIND
+
+        ! Transformation of sub-triangle center to physical coordinates
+        call ocn_init_tri_coordinate_transform(rCenter, sCenter, xTri, yTri, xCenter, yCenter)
+
+        ! Evaluate bathymetry
+        call problem_bathymetry(xCenter, yCenter, subgridBathymetryValues(nSubgridCV))
+
+        ! Transformation of sub-triangle vertices to physical coordinates
+        do i = 1,3
+           call ocn_init_tri_coordinate_transform(rSubgridPoints(i,iPt), sSubgridPoints(i,iPt), xTri, yTri, xSubgridPoints(i), ySubgridPoints(i))
+        enddo
+ 
+        ! Calculate area of sub-triangle
+        call ocn_init_tri_area(xSubgridPoints(:), ySubgridPoints(:), subgridAreas(nSubgridCV))
+
+        ! Optionally evalulate ssh !
+        if (present(subgridsshValues)) then
+           call problem_ssh(xCenter, yCenter, subgridBathymetryValues(nSubgridCV), subgridSshValues(nSubgridCV))
+        endif
+
+        if (present(subgridUValues).and.present(subgridVValues)) then
+           call problem_velocity(xCenter, yCenter, subgridUValues(nSubgridCV), subgridVValues(nSubgridCV))
+        endif   
+
+        if ( present(xSubgridCell) ) then
+            xSubgridCell(nSubgridCV) = xCenter ;
+        endif
+
+        if ( present(ySubgridCell) ) then
+            ySubgridCell(nSubgridCV) = yCenter ;
+        end if
+      enddo
+
+   !--------------------------------------------------------------------
+
+   end subroutine ocn_init_evaluate_subgrid_data!}}}
+
+!***********************************************************************
+!
+!  routine ocn_init_tri_area
+!
+!> \brief   Compute triangle area
+!> \author  Steven Brus
+!> \date    November 2022
+!> \details Compute area of triangle given vertex coordinates
+!>  
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_init_tri_area(x, y, area)!{{{
+           
+      implicit none 
+
+      real (kind=RKIND), intent(in) :: x(3), y(3)
+      real (kind=RKIND), intent(out) :: area
+
+   !--------------------------------------------------------------------
+
+      area = 0.5_RKIND*(x(2)*y(3) - y(2)*x(3) - x(1)*y(3) + y(1)*x(3) + x(1)*y(2) - y(1)*x(2))
+
+   !--------------------------------------------------------------------
+
+   end subroutine ocn_init_tri_area!}}}
+
+!***********************************************************************
+!
+!  routine ocn_init_tri_coordinate_transform
+!
+!> \brief   Transform reference triangle coordinates to mesh coordinates
+!> \author  Steven Brus
+!> \date    November 2022
+!> \details Evaluate r,s reference coordinates in x,y mesh
+!>  
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_init_tri_coordinate_transform(r, s, xTri, yTri, x, y)!{{{
+           
+      implicit none 
+
+      real (kind=RKIND), intent(in) :: r, s
+      real (kind=RKIND), intent(in) :: xTri(3), yTri(3)
+      real (kind=RKIND), intent(out) :: x, y
+
+   !--------------------------------------------------------------------
+
+      ! Transformation of sub-triangle center to physical coordinates
+      x = 0.5_RKIND*(-(r+s)*xTri(1) + (1.0_RKIND+r)*xTri(2) + (1.0_RKIND+s)*xTri(3))
+      y = 0.5_RKIND*(-(r+s)*yTri(1) + (1.0_RKIND+r)*yTri(2) + (1.0_RKIND+s)*yTri(3))
+
+   !--------------------------------------------------------------------
+
+   end subroutine ocn_init_tri_coordinate_transform!}}}
+
+!***********************************************************************
+!
+!  routine ocn_init_vertical_integration
+!
+!> \brief   Compute the wet volume per unit area lookup table
+!> \author  Steven Brus
+!> \date    November 2022
+!> \details Integrate the wet fraction over discrete ssh values 
+!>  
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_init_vertical_integration(iCV, subgridSshTableRange, nSubgridCV, subgridBathymetryValues, subgridAreas, &
+                                            subgridWetVolumeTable, subgridWetFractionTable)!{{{
+           
+      implicit none 
+
+      integer, intent(in) :: iCV
+      real (kind=RKIND), dimension(:,:), intent(in) :: subgridSshTableRange
+      integer, intent(in) :: nSubgridCV
+      real (kind=RKIND), dimension(:), intent(in) :: subgridBathymetryValues, subgridAreas
+      real (kind=RKIND), dimension(:,:), intent(inout) :: subgridWetVolumeTable
+      real (kind=RKIND), dimension(:,:), intent(inout) :: subgridWetFractionTable
+
+      real (kind=RKIND) :: deltaZ, ssh, pVal
+      integer :: lev, tri
+      
+
+   !--------------------------------------------------------------------
+   
+      deltaZ = (subgridSshTableRange(2,iCV)-subgridSshTableRange(1,iCV))/real(nSubgridTableLevels-1,RKIND);
+
+      subgridWetVolumeTable(1,iCV) = subgrid_thin_layer
+      ssh = subgridSshTableRange(1,iCV) + deltaZ
+      subgridWetVolumeTable(2,iCV) = subgridWetVolumeTable(1,iCV)*sum(subgridAreas(1:nSubgridCV))
+      subgridWetFractionTable(1,iCV) = 0.0_RKIND
+      do lev = 2,nSubgridTableLevels
+        do tri = 1,nSubgridCV
+
+          pVal = 0.0_RKIND
+          if (subgridBathymetryValues(tri) + ssh >= 0.0_RKIND) then
+            pVal = 1.0_RKIND
+          endif
+
+          subgridWetVolumeTable(lev,iCV) = subgridWetVolumeTable(lev,iCV) + pVal*deltaZ*subgridAreas(tri)
+          subgridWetFractionTable(lev,iCV) = subgridWetFractionTable(lev,iCV) + pVal*subgridAreas(tri)
+  
+        enddo
+
+        if (lev < nSubgridTableLevels) then
+          subgridWetVolumeTable(lev+1,iCV) = subgridWetVolumeTable(lev,iCV)
+        endif
+
+        subgridWetVolumeTable(lev,iCV) = subgridWetVolumeTable(lev,iCV)/sum(subgridAreas(1:nSubgridCV))
+        subgridWetFractionTable(lev,iCV) = subgridWetFractionTable(lev,iCV)/sum(subgridAreas(1:nSubgridCV))
+
+        ssh = ssh + deltaZ
+      enddo
+
+   !--------------------------------------------------------------------
+
+   end subroutine ocn_init_vertical_integration!}}}
+
+
+!***********************************************************************
+!
+!  routine ocn_init_grid_average
+!
+!> \brief   Compute thickness and velocity averages over wet area 
+!> \author  Steven Brus
+!> \date    November 2022
+!> \details Compute the thickness and velocity averages over an area
+!>          based on the subgrid wet area
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_init_grid_average(nSubgridCV, subgridBathymetryValues, subgridSshValues, subgridAreas, subgridThicknessAverage, &
+                                   subgridUValues, subgridVValues, subgridUAverage, subgridVAverage)!{{{
+           
+      implicit none 
+
+      integer, intent(in) :: nSubgridCV
+      real (kind=RKIND), dimension(:), intent(in) :: subgridBathymetryValues, subgridAreas
+      real (kind=RKIND), dimension(:), intent(in) :: subgridSshValues
+      real (kind=RKIND), intent(inout) :: subgridThicknessAverage
+      real (kind=RKIND), dimension(:), intent(in), optional :: subgridUValues, subgridVValues
+      real (kind=RKIND), intent(inout), optional :: subgridUAverage, subgridVAverage
+
+      real (kind=RKIND) :: deltaZ, ssh, pVal
+      real (kind=RKIND) :: averageDepth, cellArea, layerThicknessValue
+      real (kind=RKIND) :: H_int 
+      real (kind=RKIND) :: HU_int, HV_int
+      logical :: computeVelAverage
+      integer :: lev, tri
+      
+
+   !--------------------------------------------------------------------
+   
+      averageDepth = sum(subgridBathymetryValues(1:nSubgridCV))/real(nSubgridCV,RKIND) 
+
+      computeVelAverage = .false.
+      if (present(subgridUValues) .and. present(subgridVValues) .and. &
+          present(subgridUAverage) .and. present(subgridVAverage)) then
+         computeVelAverage = .true.
+      endif
+
+      H_int = 0.0_RKIND
+      cellArea = 0.0_RKIND
+      HU_int = 0.0_RKIND
+      HV_int = 0.0_RKIND
+      do tri = 1,nSubgridCV
+
+        layerThicknessValue = subgridBathymetryValues(tri) + subgridSshValues(tri)
+        if (layerThicknessValue <= config_drying_min_cell_height + eps) then
+          layerThicknessValue = config_drying_min_cell_height + eps
+        endif
+        H_int = H_int + layerThicknessValue*subgridAreas(tri)
+        cellArea = cellArea + subgridAreas(tri)
+
+        if (computeVelAverage) then
+           HU_int = HU_int + layerThicknessValue*subgridUValues(tri)*subgridAreas(tri)
+           HV_int = HV_int + layerThicknessValue*subgridVValues(tri)*subgridAreas(tri)
+        endif
+      enddo
+
+      if (computeVelAverage) then
+         subgridUAverage = HU_int/H_int
+         subgridVAverage = HV_int/H_int
+      endif
+      subgridThicknessAverage = H_int/cellArea
+      
+
+   !--------------------------------------------------------------------
+
+   end subroutine ocn_init_grid_average!}}}
+
+!***********************************************************************
+!
+!  routine ocn_init_wet_average_ssh
+!
+!> \brief   
+!> \author  Steven Brus
+!> \date    November 2022
+!> \details 
+!>  
+!
+!-----------------------------------------------------------------------
+   subroutine ocn_init_wet_average_ssh( nSubgridCV, subgridBathymetryValues, subgridSshValues, subgridAreas, sshWetAverage)!{{{
+      implicit none
+
+      integer, intent(in) :: nSubgridCV
+      real (kind=RKIND), dimension(:), intent(in) :: subgridBathymetryValues, subgridAreas
+      real (kind=RKIND), dimension(:), intent(in) :: subgridSshValues
+      real (kind=RKIND), intent(inout) :: sshWetAverage
+
+      real (kind=RKIND) :: deltaZ, ssh, pVal
+      real (kind=RKIND) :: averageDepth, wetArea, layerThicknessValue
+      integer :: lev, tri
+
+
+      sshWetAverage = 0.0_RKIND 
+      wetArea = 0.0_RKIND 
+
+      do tri = 1, nsubgridCV
+        layerThicknessValue = subgridBathymetryValues(tri) + subgridSshValues(tri)
+
+        if ( layerThicknessValue > subgrid_thin_layer) then 
+           sshWetAverage = sshWetAverage + subgridSshValues(tri)*subgridAreas(tri)
+           wetArea = wetArea + subgridAreas(tri)
+        endif
+      end do
+
+      if ( WetArea > 0.0_RKIND ) then
+          sshWetAverage = sshWetAverage/WetArea ;
+      else
+         sshWetAverage = -maxval(subgridBathymetryValues(1:nsubgridCV)) + subgrid_thin_layer
+      endif
+
+      return ;
+   end subroutine ocn_init_wet_average_ssh
+
+
+end module ocn_init_subgrid
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
diff --git a/components/mpas-ocean/src/ocean.cmake b/components/mpas-ocean/src/ocean.cmake
index b353c231df79..8866a8cea3d6 100644
--- a/components/mpas-ocean/src/ocean.cmake
+++ b/components/mpas-ocean/src/ocean.cmake
@@ -113,6 +113,7 @@ list(APPEND RAW_SOURCES
   core_ocean/shared/mpas_ocn_vel_tidal_potential.F
   core_ocean/shared/mpas_ocn_stokes_drift.F
   core_ocean/shared/mpas_ocn_manufactured_solution.F
+  core_ocean/shared/mpas_ocn_subgrid.F
 )
 
 set(OCEAN_DRIVER
diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile
index 5c5bc21cff12..35017551c754 100644
--- a/components/mpas-ocean/src/shared/Makefile
+++ b/components/mpas-ocean/src/shared/Makefile
@@ -69,6 +69,7 @@ OBJS = mpas_ocn_init_routines.o \
 	   mpas_ocn_time_average_coupled.o \
 	   mpas_ocn_framework_forcing.o \
 	   mpas_ocn_time_varying_forcing.o \
+	   mpas_ocn_subgrid.o \
 	   mpas_ocn_wetting_drying.o \
 	   mpas_ocn_vel_tidal_potential.o \
 	   mpas_ocn_vel_forcing_topographic_wave_drag.o \
@@ -84,7 +85,7 @@ mpas_ocn_init_routines.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o
 
 mpas_ocn_tendency.o: mpas_ocn_high_freq_thickness_hmix_del2.o mpas_ocn_tracer_surface_restoring.o mpas_ocn_thick_surface_flux.o mpas_ocn_tracer_short_wave_absorption.o mpas_ocn_tracer_advection.o mpas_ocn_tracer_hmix.o mpas_ocn_tracer_nonlocalflux.o mpas_ocn_surface_bulk_forcing.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_tracer_surface_flux_to_tend.o mpas_ocn_tracer_interior_restoring.o mpas_ocn_tracer_exponential_decay.o mpas_ocn_tracer_ideal_age.o mpas_ocn_tracer_TTD.o mpas_ocn_vmix.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_frazil_forcing.o mpas_ocn_tidal_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_tracer_CFC.o mpas_ocn_diagnostics.o mpas_ocn_wetting_drying.o mpas_ocn_vel_self_attraction_loading.o  mpas_ocn_vel_tidal_potential.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_thick_hadv.o mpas_ocn_thick_vadv.o mpas_ocn_vel_hadv_coriolis.o mpas_ocn_vel_pressure_grad.o mpas_ocn_vel_vadv.o mpas_ocn_vel_hmix.o mpas_ocn_vel_forcing.o mpas_ocn_manufactured_solution.o
 
-mpas_ocn_diagnostics.o: mpas_ocn_thick_ale.o mpas_ocn_equation_of_state.o mpas_ocn_gm.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_vertical_advection.o mpas_ocn_submesoscale_eddies.o
+mpas_ocn_diagnostics.o: mpas_ocn_thick_ale.o mpas_ocn_equation_of_state.o mpas_ocn_gm.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_vertical_advection.o mpas_ocn_submesoscale_eddies.o mpas_ocn_subgrid.o
 
 mpas_ocn_diagnostics_variables.o: mpas_ocn_config.o
 
@@ -184,7 +185,7 @@ mpas_ocn_surface_land_ice_fluxes.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_
 
 mpas_ocn_frazil_forcing.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_equation_of_state.o
 
-mpas_ocn_tidal_forcing.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_equation_of_state.o mpas_ocn_diagnostics_variables.o mpas_ocn_mesh.o
+mpas_ocn_tidal_forcing.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_equation_of_state.o mpas_ocn_diagnostics_variables.o mpas_ocn_mesh.o mpas_ocn_subgrid.o
 
 mpas_ocn_transport_tests.o: mpas_ocn_config.o
 
@@ -238,6 +239,8 @@ mpas_ocn_stokes_drift.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_diagnos
 
 mpas_ocn_manufactured_solution.o:  mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o 
 
+mpas_ocn_subgrid.o: mpas_ocn_diagnostics_variables.o mpas_ocn_config.o mpas_ocn_constants.o
+
 clean:
 	$(RM) *.o *.i *.mod *.f90
 
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F
index 77796959538d..cdb398f43ff7 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F
@@ -38,6 +38,7 @@ module ocn_diagnostics
    use ocn_mesh
    use ocn_surface_land_ice_fluxes
    use ocn_vertical_advection
+   use ocn_subgrid
 
    implicit none
    private
@@ -239,10 +240,19 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMe
 !     endif
 #endif
 
+      !
+      ! Z-coordinates
+      ! This section must be placed in the code before computing the density.
+      !
+      ! inputs : layerThickness
+      ! outputs : zMid, zTop, ssh)
+      call ocn_diagnostic_solve_z_coordinates(layerThickness, zMid, zTop, ssh)
+
       ! inputs: layerThickness, normalVelocity
       ! output: layerThickEdgeMean, layerThickEdgeDrag, layerThickEdgeFlux
       call ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, &
-                                                   layerThickness, restingThickness)
+                                                   layerThickness, restingThickness, &
+                                                   ssh)
 
       ! inputs: normalVelocity
       ! outputs: relativeVorticity, circulation
@@ -337,13 +347,6 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMe
       call ocn_diagnostic_solve_surface_pressure(forcingPool, atmosphericPressure, &
                 seaIcePressure, surfacePressure)
 
-      !
-      ! Z-coordinates
-      ! This section must be placed in the code before computing the density.
-      !
-      ! inputs : layerThickness
-      ! outputs : zMid, zTop, ssh)
-      call ocn_diagnostic_solve_z_coordinates(layerThickness, zMid, zTop, ssh)
 
       !
       ! equation of state
@@ -807,7 +810,8 @@ end subroutine ocn_diagnostic_solve_wctEdge!}}}
 !-----------------------------------------------------------------------
 
    subroutine ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, &
-                                                      layerThickness, restingThickness)!{{{
+                                                      layerThickness, restingThickness, &
+                                                      ssh)!{{{
 
       !-----------------------------------------------------------------
       ! input variables
@@ -817,6 +821,8 @@ subroutine ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, &
          normalVelocity,  &!< [in] transport
          layerThickness,  &!< [in] layer thickness at cell center
          restingThickness  !< [in] initial layer thickness at cell center
+      real (kind=RKIND), dimension(:), intent(in) :: &
+         ssh               !< [in] sea surface height at cell center
 
       !-----------------------------------------------------------------
       ! output variables
@@ -879,64 +885,79 @@ subroutine ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, &
       case (thickEdgeFluxCenter)
          ! Use centered (mean) thickness as flux value
 
+         if ( config_use_subgrid_wetting_drying ) then
+
+            call ocn_subgrid_layerThickEdgeFlux_center(ssh, layerThickEdgeFlux)
+
+         else
+
 #ifdef MPAS_OPENACC
-         !$acc parallel loop collapse(2) &
-         !$acc    present(layerThickEdgeFlux, layerThickEdgeMean)
+            !$acc parallel loop collapse(2) &
+            !$acc    present(layerThickEdgeFlux, layerThickEdgeMean)
 #else
-         !$omp parallel
-         !$omp do schedule(runtime) private(k)
+            !$omp parallel
+            !$omp do schedule(runtime) private(k)
 #endif
-         do iEdge = 1, nEdgesAll
-         do k = 1,nVertLevels
-            layerThickEdgeFlux(k,iEdge) = &
-            layerThickEdgeMean(k,iEdge)
-         end do
-         end do
+            do iEdge = 1, nEdgesAll
+            do k = 1,nVertLevels
+               layerThickEdgeFlux(k,iEdge) = &
+               layerThickEdgeMean(k,iEdge)
+            end do
+            end do
 #ifndef MPAS_OPENACC
-         !$omp end do
-         !$omp end parallel
+            !$omp end do
+            !$omp end parallel
 #endif
 
+         end if  
+
       case (thickEdgeFluxUpwind)
          ! Use upwind thickness as the edge flux value
 
+         if (config_use_subgrid_wetting_drying) then
+
+            call ocn_subgrid_layerThickEdgeFlux_upwind(ssh, normalVelocity, layerThickness, layerThickEdgeFlux)
+
+         else
 #ifdef MPAS_OPENACC
-         !$acc parallel loop &
-         !$acc    present(normalVelocity, layerThickness, &
-         !$acc            minLevelEdgeBot, maxLevelEdgeTop, &
-         !$acc            layerThickEdgeFlux, cellsOnEdge) &
-         !$acc    private(k, kmin, kmax, cell1, cell2)
+            !$acc parallel loop &
+            !$acc    present(normalVelocity, layerThickness, &
+            !$acc            minLevelEdgeBot, maxLevelEdgeTop, &
+            !$acc            layerThickEdgeFlux, cellsOnEdge) &
+            !$acc    private(k, kmin, kmax, cell1, cell2)
 #else
-         !$omp parallel
-         !$omp do schedule(runtime) private(k, kmin, kmax, cell1, cell2)
+            !$omp parallel
+            !$omp do schedule(runtime) private(k, kmin, kmax, cell1, cell2)
 #endif
-         do iEdge = 1, nEdgesAll
-            kmin = minLevelEdgeBot(iEdge)
-            kmax = maxLevelEdgeTop(iEdge)
-            cell1 = cellsOnEdge(1,iEdge)
-            cell2 = cellsOnEdge(2,iEdge)
-            do k=1,nVertLevels
-               ! initialize layerThicknessEdgeFlux to avoid divide by
-               ! zero and NaN problems.
-               layerThickEdgeFlux(k,iEdge) = -1.0e34_RKIND
-            end do
-            do k = kmin,kmax
-               if (normalVelocity(k,iEdge) > 0.0_RKIND) then
-                  layerThickEdgeFlux(k,iEdge) = layerThickness(k,cell1)
-               elseif (normalVelocity(k,iEdge) < 0.0_RKIND) then
-                  layerThickEdgeFlux(k,iEdge) = layerThickness(k,cell2)
-               else
-                  layerThickEdgeFlux(k,iEdge) = &
-                                      max(layerThickness(k,cell1), &
-                                          layerThickness(k,cell2))
-               end if
+            do iEdge = 1, nEdgesAll
+               kmin = minLevelEdgeBot(iEdge)
+               kmax = maxLevelEdgeTop(iEdge)
+               cell1 = cellsOnEdge(1,iEdge)
+               cell2 = cellsOnEdge(2,iEdge)
+               do k=1,nVertLevels
+                  ! initialize layerThicknessEdgeFlux to avoid divide by
+                  ! zero and NaN problems.
+                  layerThickEdgeFlux(k,iEdge) = -1.0e34_RKIND
+               end do
+               do k = kmin,kmax
+                  if (normalVelocity(k,iEdge) > 0.0_RKIND) then
+                     layerThickEdgeFlux(k,iEdge) = layerThickness(k,cell1)
+                  elseif (normalVelocity(k,iEdge) < 0.0_RKIND) then
+                     layerThickEdgeFlux(k,iEdge) = layerThickness(k,cell2)
+                  else
+                     layerThickEdgeFlux(k,iEdge) = &
+                                         max(layerThickness(k,cell1), &
+                                             layerThickness(k,cell2))
+                  end if
+               end do
             end do
-         end do
 #ifndef MPAS_OPENACC
-         !$omp end do
-         !$omp end parallel
+            !$omp end do
+            !$omp end parallel
 #endif
 
+        end if
+
       case (thickEdgeFluxConstant)
          ! Use linearized version H*u where H is constant in time
 
@@ -2600,8 +2621,21 @@ subroutine ocn_diagnostic_solve_z_coordinates(layerThickness, zMid, zTop, ssh)!{
               + layerThickness(k  ,iCell)
          end do
 
-         ! copy zTop(1,iCell) into sea-surface height array
-         ssh(iCell) = zTop(minLevelCell(iCell),iCell)
+         if (config_use_subgrid_wetting_drying) then
+            call ocn_subgrid_ssh_lookup(layerThickness(1,iCell), &
+                                        subgridWetVolumeCellTable(:,iCell), &
+                                        subgridSshCellTableRange(:,iCell), &
+                                        bottomDepth(iCell), &
+                                        subgridCellBathymetryMin(iCell), & 
+                                        ssh(iCell))
+
+
+            zTop(1,iCell) = ssh(iCell)
+            zMid(1,iCell) = -bottomDepth(iCell) + 0.5_RKIND*layerThickness(1,iCell)         
+         else
+            ! copy zTop(1,iCell) into sea-surface height array
+            ssh(iCell) = zTop(minLevelCell(iCell),iCell)
+         end if
 
       end do
 #ifndef MPAS_OPENACC
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
index 567c3efc183c..2b275a6f0811 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
@@ -140,6 +140,9 @@ module ocn_diagnostics_variables
    real (kind=RKIND), dimension(:,:), pointer :: wettingVelocityFactor
    type (field2DReal), pointer :: wettingVelocityField
 
+   real (kind=RKIND), dimension(:,:), pointer :: layerThicknessCellWetDry
+   real (kind=RKIND), dimension(:), pointer :: sshCellWetDry
+
    real (kind=RKIND), dimension(:,:), pointer :: vertDiffTopOfCell
    real (kind=RKIND), dimension(:,:,:), pointer :: vertNonLocalFlux
    integer, dimension(:,:), pointer :: rediLimiterCount
@@ -262,6 +265,10 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
                    wettingVelocityFactor)
          call mpas_pool_get_field(diagnosticsPool, 'wettingVelocityFactor', &
                    wettingVelocityField)
+         call mpas_pool_get_array(diagnosticsPool, 'layerThicknessCellWetDry', &
+                   layerThicknessCellWetDry)
+         call mpas_pool_get_array(diagnosticsPool, 'sshCellWetDry', &
+                   sshCellWetDry)
       end if
 
       if (config_use_GM.or.config_use_Redi) then
@@ -1246,6 +1253,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
       end if
       if ( trim(config_ocean_run_mode) == 'forward' ) then
          nullify(wettingVelocityFactor)
+         nullify(layerThicknessCellWetDry)
+         nullify(sshCellWetDry)
       end if
       if (config_use_GM.or.config_use_Redi) then
          nullify(RediKappa,         &
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F b/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F
new file mode 100644
index 000000000000..82208a971687
--- /dev/null
+++ b/components/mpas-ocean/src/shared/mpas_ocn_subgrid.F
@@ -0,0 +1,589 @@
+! Copyright (c) 2013,  Los Alamos National Security, LLC (LANS)
+! and the University Corporation for Atmospheric Research (UCAR).
+!
+! Unless noted otherwise source code is licensed under the BSD license.
+! Additional copyright and license information can be found in the LICENSE file
+! distributed with this code, or at http://mpas-dev.github.com/license.html
+!
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  ocn_subgrid
+!
+!> \brief   MPAS ocean subgrid wetting and drying
+!> \authors Steven Brus, Damrongsak Wirasaet
+!> \date    September 2022
+!> \details
+!>  This module contains routines for subgrid wetting and drying
+!>  These routines primarily determine the relationship between the
+!>  wet volume per unit area (which can be thought of a as a layer thickness
+!>  for a completely wet cell) of a control volume and a the ssh
+!>  defined over the wet fraction of that control volume. When a
+!>  control volume is wet, this relationship is linear, i.e. H = \eta + b.
+!>  However, for partially wet cells, the nonlinear relationship between
+!>  these values can be pre-computed based on subgrid scale high-resolution data.
+!>  The routines in this module are responsible for providing interpolated
+!>  values for the forward and inverse relationship to be used in
+!>  computing the tendency terms from pre-computed look-up tables.
+!>
+!>  See other details in Kennedy et al. (2019):
+!>  https://doi.org/10.1016/j.ocemod.2019.101491
+!
+!-----------------------------------------------------------------------
+
+module ocn_subgrid
+
+   use mpas_kind_types
+   use mpas_constants
+   use mpas_derived_types
+   use mpas_pool_routines
+   use mpas_timekeeping
+   use mpas_timer
+   use ocn_constants
+   use ocn_config
+   use ocn_mesh
+   use ocn_diagnostics_variables
+
+   implicit none
+   private
+   save
+
+   !--------------------------------------------------------------------
+   !
+   ! Public parameters
+   !
+   !--------------------------------------------------------------------
+
+   real(kind=RKIND), public, dimension(:,:), pointer :: &
+      subgridWetVolumeCellTable, &
+      subgridWetVolumeEdgeTable, &
+      subgridWetVolumeVertexTable, &
+      subgridWetFractionCellTable, &
+      subgridWetFractionEdgeTable, &
+      subgridWetFractionVertexTable, &
+      subgridSshCellTableRange, &
+      subgridSshEdgeTableRange, &
+      subgridSshVertexTableRange
+   real(kind=RKIND), public, dimension(:), pointer :: &
+      subgridEdgeBathymetryMean, &
+      subgridVertexBathymetryMean, &
+      subgridCellBathymetryMin, &
+      subgridEdgeBathymetryMin, &
+      subgridVertexBathymetryMin, &
+      subgridLayerThicknessDebug
+
+   integer, public, pointer :: nSubgridTableLevels
+
+   !--------------------------------------------------------------------
+   !
+   ! Public member functions
+   !
+   !--------------------------------------------------------------------
+
+   public :: ocn_subgrid_layerThickEdgeFlux_center, &
+             ocn_subgrid_layerThickEdgeFlux_upwind, &
+             ocn_subgrid_vorticity
+
+   public :: ocn_subgrid_layer_thickness_lookup, &
+             ocn_subgrid_wet_fraction_lookup, &
+             ocn_subgrid_ssh_lookup, &
+             ocn_subgrid_init
+
+   !--------------------------------------------------------------------
+   !
+   ! Private module variables
+   !
+   !--------------------------------------------------------------------
+
+
+
+contains
+
+!***********************************************************************
+!
+!  routine ocn_subgrid_layerThickEdgeFlux_center
+!
+!> \brief   Copmutes centered layer thickness edge flux using subgrid info
+!> \author  Steven Brus, Damrongsak Wirasaet
+!> \date    September 2022
+!> \details
+!>
+!>
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_subgrid_layerThickEdgeFlux_center(ssh, layerThickEdgeFlux)!{{{
+
+      !-----------------------------------------------------------------
+      ! input variables
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:), intent(in) :: ssh
+
+      !-----------------------------------------------------------------
+      ! output variables
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:,:), intent(out) :: layerThickEdgeFlux
+
+      !-----------------------------------------------------------------
+      ! local variables
+      !-----------------------------------------------------------------
+
+      integer :: iEdge
+      integer :: cell1, cell2
+      real (kind=RKIND) :: sshMean
+      real(kind=RKIND):: eps = 1.0e-10_RKIND
+
+
+      do iEdge = 1, nEdgesAll
+         cell1 = cellsOnEdge(1,iEdge)
+         cell2 = cellsOnEdge(2,iEdge)
+         sshMean = 0.5_RKIND * (ssh(cell1) + ssh(cell2))
+         call ocn_subgrid_layer_thickness_lookup(sshMean, &
+                                                 subgridWetVolumeEdgeTable(:,iEdge), &
+                                                 subgridSshEdgeTableRange(:,iEdge), &
+                                                 subgridEdgeBathymetryMean(iEdge), &
+                                                 layerThickEdgeFlux(1,iEdge) )
+
+         if ( layerThickEdgeFlux(1,iEdge) < eps ) then
+             layerThickEdgeFlux(1,iEdge) = layerThickEdgeMean(1,iEdge) ;
+         end if
+      end do
+
+   end subroutine ocn_subgrid_layerThickEdgeFlux_center!}}}
+
+!***********************************************************************
+!
+!  routine ocn_subgrid_layerThickEdgeFlux_upwind
+!
+!> \brief   Copmutes upwind layer thickness edge flux using subgrid info
+!> \author  Steven Brus, Damrongsak Wirasaet
+!> \date    September 2022
+!> \details
+!>
+!>
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_subgrid_layerThickEdgeFlux_upwind(ssh, normalVelocity, layerThickness, layerThickEdgeFlux)!{{{
+
+      !-----------------------------------------------------------------
+      ! input variables
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:), intent(in) :: ssh
+      real (kind=RKIND), dimension(:,:), intent(in) :: normalVelocity
+      real (kind=RKIND), dimension(:,:), intent(in) :: layerThickness
+
+      !-----------------------------------------------------------------
+      ! output variables
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:,:), intent(out) :: layerThickEdgeFlux
+
+      !-----------------------------------------------------------------
+      ! local variables
+      !-----------------------------------------------------------------
+
+      integer :: iEdge
+      integer :: cell1, cell2
+      real (kind=RKIND) :: sshMean
+      real(kind=RKIND):: eps = 1.0e-10_RKIND
+
+      do iEdge = 1, nEdgesAll
+         cell1 = cellsOnEdge(1,iEdge)
+         cell2 = cellsOnEdge(2,iEdge)
+         if (normalVelocity(1,iEdge) > 0.0_RKIND) then
+            call ocn_subgrid_layer_thickness_lookup(ssh(cell1), &
+                                                    subgridWetVolumeEdgeTable(:,iEdge), &
+                                                    subgridSshEdgeTableRange(:,iEdge), &
+                                                    subgridEdgeBathymetryMean(iEdge), &
+                                                    layerThickEdgeFlux(1,iEdge))
+         elseif (normalVelocity(1,iEdge) < 0.0_RKIND) then
+            call ocn_subgrid_layer_thickness_lookup(ssh(cell2), &
+                                                    subgridWetVolumeEdgeTable(:,iEdge), &
+                                                    subgridSshEdgeTableRange(:,iEdge), &
+                                                    subgridEdgeBathymetryMean(iEdge), &
+                                                    layerThickEdgeFlux(1,iEdge))
+         else
+            sshMean = 0.5_RKIND*(ssh(cell1) + ssh(cell2))
+            call ocn_subgrid_layer_thickness_lookup(sshMean, &
+                                                    subgridWetVolumeEdgeTable(:,iEdge), &
+                                                    subgridSshEdgeTableRange(:,iEdge), &
+                                                    subgridEdgeBathymetryMean(iEdge), &
+                                                    layerThickEdgeFlux(1,iEdge))
+         end if
+
+         if ( layerThickEdgeFlux(1,iEdge) < eps ) then
+            cell1 = cellsOnEdge(1,iEdge)
+            cell2 = cellsOnEdge(2,iEdge)
+
+            layerThickEdgeFlux(1,iEdge) = 0.5_RKIND * &
+                                     ( layerThickness(1,cell1) + &
+                                       layerThickness(1,cell2) )
+         end if
+      end do
+
+   end subroutine ocn_subgrid_layerThickEdgeFlux_upwind!}}}
+
+!***********************************************************************
+!
+!  routine ocn_subgrid_vorticity
+!
+!> \brief   Copmutes vorticity using subgrid info
+!> \author  Steven Brus, Damrongsak Wirasaet
+!> \date    September 2022
+!> \details
+!>
+!>
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_subgrid_vorticity(ssh, normalizedRelativeVorticityVertex, normalizedPlanetaryVorticityVertex)!{{{
+
+      !-----------------------------------------------------------------
+      ! input variables
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:), intent(in) :: ssh
+
+      !-----------------------------------------------------------------
+      ! output variables
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:,:), intent(out) :: normalizedRelativeVorticityVertex
+      real (kind=RKIND), dimension(:,:), intent(out) :: normalizedPlanetaryVorticityVertex
+
+      !-----------------------------------------------------------------
+      ! local variables
+      !-----------------------------------------------------------------
+
+      integer :: iVertex, vertex
+      integer :: i, k
+      real (kind=RKIND) :: sshVertex
+      real (kind=RKIND) :: invAreaTri1
+      real (kind=RKIND) :: layerThicknessVertex
+
+      do iVertex = 1, nVerticesHalo(2)
+         invAreaTri1 = 1.0_RKIND / areaTriangle(iVertex)
+         do k = 1, maxLevelVertexBot(iVertex)
+            sshVertex = 0.0_RKIND
+            do i = 1, vertexDegree
+               sshVertex = sshVertex + ssh(cellsOnVertex(i,iVertex)) &
+                                    * kiteAreasOnVertex(i,iVertex)
+            end do
+            sshVertex = sshVertex * invAreaTri1
+
+            call ocn_subgrid_layer_thickness_lookup(sshVertex, &
+                                                    subgridWetVolumeVertexTable(:,iVertex), &
+                                                    subgridSshVertexTableRange(:,iVertex), &
+                                                    subgridVertexBathymetryMean(iVertex), &
+                                                    layerThicknessVertex)
+            if (layerThicknessVertex == 0) cycle
+
+            normalizedRelativeVorticityVertex(k,iVertex) = relativeVorticity(k,iVertex) / layerThicknessVertex
+            normalizedPlanetaryVorticityVertex(k,iVertex) = fVertex(iVertex) / layerThicknessVertex
+         end do
+      end do
+
+   end subroutine ocn_subgrid_vorticity!}}}
+
+!***********************************************************************
+!
+!  routine ocn_subgrid_layer_thickness_lookup
+!
+!> \brief   Forward subgrid lookup (ssh -> layerThicknes)
+!> \author  Steven Brus, Damrongsak Wirasaet
+!> \date    September 2022
+!> \details Returns the wet volume per unit area (layerThick) for a
+!>          given ssh (defined over the wet-fraction based on subgrid
+!>          lookup table information
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_subgrid_layer_thickness_lookup(zeta, &
+                                                 subgridTable, &
+                                                 subgridTableRange, &
+                                                 bathymetry, &
+                                                 layerThick)!{{{
+
+      !-----------------------------------------------------------------
+      ! input variables
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), intent(in) :: zeta
+      real (kind=RKIND), dimension(:), intent(in) :: subgridTable
+      real (kind=RKIND), dimension(:), intent(in) :: subgridTableRange
+      real (kind=RKIND), intent(in) :: bathymetry
+
+      !-----------------------------------------------------------------
+      ! output variables
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), intent(out) :: layerThick
+
+      !-----------------------------------------------------------------
+      ! local variables
+      !-----------------------------------------------------------------
+
+      integer :: lev
+
+      real (kind=RKIND) :: tableMin
+      real (kind=RKIND) :: tableMax
+      real (kind=RKIND) :: deltaZ
+      real (kind=RKIND) :: zeta0, zeta1
+
+      tableMin = subgridTableRange(1)
+      tableMax = subgridTableRange(2)
+      deltaZ = (tableMax - tableMin)/real(nSubgridTableLevels-1,RKIND)
+
+      if (zeta >= tableMax) then
+        layerThick = zeta + bathymetry
+      else if (zeta <= tableMin) then
+        layerThick = 0.0_RKIND
+      else
+        do lev = 1, nSubgridTableLevels-1
+          zeta0 = (real(lev,RKIND)-1.0_RKIND)*deltaZ + tableMin
+          zeta1 = zeta0 + deltaZ
+
+          if ((zeta <= zeta1) .and. (zeta >= zeta0)) then
+            layerThick = ((zeta-zeta0)*subgridTable(lev+1) - (zeta-zeta1)*subgridTable(lev))/deltaZ
+            return
+          end if
+
+        end do
+      end if
+
+   end subroutine ocn_subgrid_layer_thickness_lookup!}}}
+
+!***********************************************************************
+!
+!  routine ocn_subgrid_wet_fraction_lookup
+!
+!> \brief   Wet fraction lookup
+!> \author  Steven Brus, Damrongsak Wirasaet
+!> \date    September 2022
+!> \details Returns the wet fraction for a given ssh (defined over the
+!>          wet fraction) value based on subgrid lookup table information
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_subgrid_wet_fraction_lookup(zeta, &
+                                              subgridTable, &
+                                              subgridTableRange, &
+                                              wetFraction)!{{{
+
+      !-----------------------------------------------------------------
+      ! input variables
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), intent(in) :: zeta
+      real (kind=RKIND), dimension(:), intent(in) :: subgridTable
+      real (kind=RKIND), dimension(:), intent(in) :: subgridTableRange
+
+      !-----------------------------------------------------------------
+      ! output variables
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), intent(out) :: wetFraction
+
+      !-----------------------------------------------------------------
+      ! local variables
+      !-----------------------------------------------------------------
+
+      integer :: lev
+
+      real (kind=RKIND) :: tableMin
+      real (kind=RKIND) :: tableMax
+      real (kind=RKIND) :: deltaZ
+      real (kind=RKIND) :: zeta0, zeta1
+
+      tableMin = subgridTableRange(1)
+      tableMax = subgridTableRange(2)
+      deltaZ = (tableMax - tableMin)/real(nSubgridTableLevels-1,RKIND)
+
+      if (zeta >= tableMax) then
+        wetFraction = 1.0_RKIND
+      else if (zeta <= tableMin) then
+        wetFraction = 0.0_RKIND
+      else
+        do lev = 1, nSubgridTableLevels-1
+          zeta0 = (real(lev,RKIND)-1.0_RKIND)*deltaZ + tableMin
+          zeta1 = zeta0 + deltaZ
+
+          if ((zeta <= zeta1) .and. (zeta >= zeta0)) then
+            wetFraction = ((zeta-zeta0)*subgridTable(lev+1) - (zeta-zeta1)*subgridTable(lev))/deltaZ
+            return
+          end if
+
+        end do
+      end if
+
+   end subroutine ocn_subgrid_wet_fraction_lookup!}}}
+
+!***********************************************************************
+!
+!  routine ocn_subgrid_ssh_lookup
+!> \brief   Inverse subgrid lookup (layerThickness -> ssh)
+!> \author  Steven Brus, Damrongsak Wirasaet
+!> \date    September 2022
+!> \details Returns the ssh (defined over the wet fraction) for a given
+!>          wet volume per unit area value (layerThick) based on subgrid
+!>          lookup table information
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_subgrid_ssh_lookup(layerThick, &
+                                     subgridTable, &
+                                     subgridTableRange, &
+                                     bathymetryMean, &
+                                     bathymetryMin, &
+                                     zeta)!{{{
+
+      !-----------------------------------------------------------------
+      ! input variables
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), intent(in) :: layerThick
+      real (kind=RKIND), dimension(:), intent(in) :: subgridTable
+      real (kind=RKIND), dimension(:), intent(in) :: subgridTableRange
+      real (kind=RKIND), intent(in) :: bathymetryMean
+      real (kind=RKIND), intent(in) :: bathymetryMin
+
+      !-----------------------------------------------------------------
+      ! output variables
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), intent(inout) :: zeta
+
+      !-----------------------------------------------------------------
+      ! local variables
+      !-----------------------------------------------------------------
+
+      integer :: lev
+
+      real (kind=RKIND) :: deltaZ
+      real (kind=RKIND) :: tableMin, tableMax
+      real (kind=RKIND) :: layerThickMin, layerThickMax
+      real (kind=RKIND) :: layerThick0, layerThick1
+      real (kind=RKIND) :: zeta0, zeta1
+      real (kind=RKIND) :: phi0, phi1
+
+      tableMin = subgridTableRange(1)
+      tableMax = subgridTableRange(2)
+      deltaZ = (tableMax - tableMin)/real(nSubgridTableLevels-1,RKIND)
+
+      layerThickMin = subgridTable(1)
+      layerThickMax = subgridTable(nSubgridTableLevels)
+
+      if (layerThick >= layerThickMax) then
+        zeta = layerThick - bathymetryMean
+
+      else if (layerThick <= layerThickMin) then
+        zeta = - bathymetryMin ! prevent_drying likely fails to ensure positive water columbn
+      else
+
+       do lev = 1, nSubgridTableLevels-1
+          zeta0 = (real(lev,RKIND)-1.0_RKIND)*deltaZ + tableMin
+          zeta1 = zeta0 + deltaZ
+
+          layerThick0 = subgridTable(lev)
+          layerThick1 = subgridTable(lev+1)
+
+          if ((layerThick <= layerThick1) .and. (layerThick >= layerThick0)) then
+            phi0 = (layerThick-layerThick1)/(layerThick0-layerThick1)
+            phi1 = (layerThick-layerThick0)/(layerThick1-layerThick0)
+            zeta = phi0*zeta0 + phi1*zeta1
+            return
+          end if
+
+        end do
+      end if
+
+   end subroutine ocn_subgrid_ssh_lookup!}}}
+
+!***********************************************************************
+!
+!  routine ocn_subgrid_init
+!
+!> \brief   Initializes subgrid wetting and drying module.
+!> \author  Steven Brus, Damrongsak Wirasaet
+!> \date    September 2022
+!> \details
+!>  This routine initializes the subgrid wetting and drying module
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_subgrid_init(domain,err)!{{{
+
+      !-----------------------------------------------------------------
+      ! input/output variables
+      !-----------------------------------------------------------------
+
+      type (domain_type), intent(inout) :: domain
+
+      !-----------------------------------------------------------------
+      ! output variables
+      !-----------------------------------------------------------------
+
+      integer, intent(out) :: err
+
+      !-----------------------------------------------------------------
+      ! local variables
+      !-----------------------------------------------------------------
+
+      type(block_type), pointer :: block
+      type(mpas_pool_type), pointer :: meshPool
+
+      err = 0
+      if (.not. config_use_subgrid_wetting_drying) then
+        return
+      end if
+
+      if ((config_ocean_run_mode == 'forward') .and. (nVertLevels .ne. 1)) then
+        call mpas_log_write('config_config_use_subgrid_wetting_drying = .true. requires single layer' , MPAS_LOG_CRIT)
+      end if
+
+      block => domain%blocklist
+      call mpas_pool_get_subpool(block%structs, 'mesh', meshPool)
+
+      call mpas_pool_get_array(meshPool, 'subgridWetVolumeCellTable', &
+                               subgridWetVolumeCellTable)
+      call mpas_pool_get_array(meshPool, 'subgridWetVolumeEdgeTable', &
+                               subgridWetVolumeEdgeTable)
+      call mpas_pool_get_array(meshPool, 'subgridWetVolumeVertexTable', &
+                               subgridWetVolumeVertexTable)
+      call mpas_pool_get_array(meshPool, 'subgridWetFractionCellTable', &
+                               subgridWetFractionCellTable)
+      call mpas_pool_get_array(meshPool, 'subgridWetFractionEdgeTable', &
+                               subgridWetFractionEdgeTable)
+      call mpas_pool_get_array(meshPool, 'subgridWetFractionVertexTable', &
+                               subgridWetFractionVertexTable)
+      call mpas_pool_get_array(meshPool, 'subgridSshCellTableRange', &
+                               subgridSshCellTableRange)
+      call mpas_pool_get_array(meshPool, 'subgridSshEdgeTableRange', &
+                               subgridSshEdgeTableRange)
+      call mpas_pool_get_array(meshPool, 'subgridSshVertexTableRange', &
+                               subgridSshVertexTableRange)
+      call mpas_pool_get_array(meshPool, 'subgridEdgeBathymetryMean', &
+                               subgridEdgeBathymetryMean)
+      call mpas_pool_get_array(meshPool, 'subgridVertexBathymetryMean', &
+                               subgridVertexBathymetryMean)
+      call mpas_pool_get_array(meshPool, 'subgridCellBathymetryMin', &
+                               subgridCellBathymetryMin)
+      call mpas_pool_get_array(meshPool, 'subgridEdgeBathymetryMin', &
+                               subgridEdgeBathymetryMin)
+      call mpas_pool_get_array(meshPool, 'subgridVertexBathymetryMin', &
+                               subgridVertexBathymetryMin)
+      call mpas_pool_get_array(meshPool, 'subgridLayerThicknessDebug', &
+                               subgridLayerThicknessDebug)
+
+      call mpas_pool_get_dimension(meshPool, 'nSubgridTableLevels', &
+                                     nSubgridTableLevels)
+
+   end subroutine ocn_subgrid_init!}}}
+
+end module ocn_subgrid
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F b/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F
index e2c1705c39b8..2a9b400dd90b 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_tidal_forcing.F
@@ -29,6 +29,7 @@ module ocn_tidal_forcing
    use ocn_config
    use ocn_mesh
    use ocn_diagnostics_variables
+   use ocn_subgrid
 
    implicit none
    private
@@ -260,7 +261,11 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo
         if (config_use_wetting_drying .and. tidalInputMask(iCell) == 1.0_RKIND) then
           ! ensure that tidal height can't force below total minimum thickness
           ! condition wrong to ensure that there isn't any drying according to criteria
-          tidalHeight = max(-bottomDepth(iCell) + (float(maxLevelCell(iCell))+1.0_RKIND)*config_drying_min_cell_height, tidalHeight)
+          if ( config_use_subgrid_wetting_drying ) then
+            tidalHeight = max(tidalHeight, subgridSShCellTableRange(3,iCell) ) ;      
+          else
+           tidalHeight = max(-bottomDepth(iCell) + (float(maxLevelCell(iCell))+1.0_RKIND)*config_drying_min_cell_height, tidalHeight)
+          end if
         end if
 
         ! compute total depth for relative thickness contribution
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F
index 8e23913960f5..b8f83cd01dce 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F
@@ -39,6 +39,19 @@ module ocn_vel_pressure_grad
    !
    !--------------------------------------------------------------------
 
+   integer, public :: &
+      pGradType   ! id for pressure gradient method selected
+
+   integer, parameter, public :: &! ids for supported methods
+      pGradTypeNone        = 0, &! none selected
+      pGradTypeSSHgrad     = 1, &! ssh gradient
+      pGradTypePZmid       = 2, &! pressure and zMid
+      pGradTypeMontPot     = 3, &! Montgomery potential
+      pGradTypeMontPotDens = 4, &! Montgomery potential and density
+      pGradTypeJacobDens   = 5, &! Jacobian from density
+      pGradTypeJacobTS     = 6, &! Jacobian from T,S
+      pGradTypeConstForced = 7   ! constant forced
+
    !--------------------------------------------------------------------
    !
    ! Public member functions
@@ -58,19 +71,6 @@ module ocn_vel_pressure_grad
       pGradOff,         &! flag for turning pressure gradient on/off
       timeIntegratorLTS  ! flag for Local Time Stepping
 
-   integer :: &
-      pGradType   ! id for pressure gradient method selected
-
-   integer, parameter ::        &! ids for supported methods
-      pGradTypeNone        = 0, &! none selected
-      pGradTypeSSHgrad     = 1, &! ssh gradient
-      pGradTypePZmid       = 2, &! pressure and zMid
-      pGradTypeMontPot     = 3, &! Montgomery potential
-      pGradTypeMontPotDens = 4, &! Montgomery potential and density
-      pGradTypeJacobDens   = 5, &! Jacobian from density
-      pGradTypeJacobTS     = 6, &! Jacobian from T,S
-      pGradTypeConstForced = 7   ! constant forced
-
    real (kind=RKIND) :: &! precomputed constants for efficiency
       density0Inv,      &! 1/density0
       gdensity0Inv,     &! g/density0
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F
index 92cb65b2bbaf..16d052df0288 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F
@@ -33,6 +33,8 @@ module ocn_wetting_drying
    use ocn_diagnostics_variables
    use ocn_gm
    use ocn_mesh
+   use ocn_subgrid
+   use ocn_vel_pressure_grad
 
    implicit none
    private
@@ -200,7 +202,7 @@ end subroutine ocn_wetting_drying_verify !}}}
 !
 !-----------------------------------------------------------------------
 
-   subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying_velocity, err) !{{{
+   subroutine ocn_prevent_drying_rk4(domain, block, dt, rkSubstepWeight, config_zero_drying_velocity, err) !{{{
 
      !-----------------------------------------------------------------
      !
@@ -208,6 +210,7 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying
      !
      !-----------------------------------------------------------------
 
+     type (domain_type), intent(inout) :: domain
      type (block_type), intent(in) :: block
      real (kind=RKIND), intent(in) :: dt
      real (kind=RKIND), intent(in) :: rkSubstepWeight
@@ -263,10 +266,10 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying
      !$omp end do
      !$omp end parallel
 
-     ! ensure cells stay wet by selectively damping cells with a damping tendency to make 
+     ! ensure cells stay wet by selectively damping cells with a damping tendency to make
      ! sure tendency doesn't dry cells
 
-     call ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
+     call ocn_wetting_drying_wettingVelocity(domain, layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
                                              normalTransportVelocity, rkSubstepWeight, wettingVelocityFactor, err)
 
      ! prevent drying from happening with selective wettingVelocityFactor
@@ -304,9 +307,9 @@ end subroutine ocn_prevent_drying_rk4 !}}}
 !>  to prevent cells from drying.
 !
 !-----------------------------------------------------------------------
-   subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
 
-       normalVelocity, dt, wettingVelocityFactor, err)!{{{
+   subroutine ocn_wetting_drying_wettingVelocity(domain, layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
+                                                 normalVelocity, dt, wettingVelocityFactor, err)!{{{
 
       !-----------------------------------------------------------------
       !
@@ -335,6 +338,10 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness
       !
       !-----------------------------------------------------------------
 
+      type (domain_type), optional, intent(inout) :: &
+         domain !< Input/Output: domain information, needed for halo exchange
+
+
       real (kind=RKIND), dimension(:,:), intent(inout) :: &
          wettingVelocityFactor          !< Input/Output: velocity wettingVelocityFactor
 
@@ -355,21 +362,39 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness
       integer :: cell1, cell2, iEdge, iCell, k, i
 
       real (kind=RKIND) :: divOutFlux
-      real (kind=RKIND) :: layerThickness
+      real (kind=RKIND) :: columnThickness
+      real (kind=RKIND) :: hCrit, hRampMin, hRampMax, hEdgeTotal
 
       character (len=100) :: log_string
 
+      integer:: cellDummy(2), cellCur, CellNei
+      real (kind=RKIND) :: sshCur, sshNei
+      real (kind=RKIND), dimension(:), pointer :: sshCell
+      real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCell
+
+
       err = 0
 
       if (.not. config_zero_drying_velocity) return
 
+      if (config_use_ssh_gradient_wetting_drying .and. &
+          pGradType /= pGradTypeSSHgrad) then
+
+         call mpas_log_write("config_use_ssh_gradient_wetting_drying requires " // &
+                             "config_pressure_gradient_type = 'ssh_gradient'", MPAS_LOG_CRIT)
+      endif
+
+      hRampMin = config_zero_drying_velocity_ramp_hmin
+      hRampMax = config_zero_drying_velocity_ramp_hmax
+      hCrit = config_drying_min_cell_height + config_drying_safety_height
+
+      layerThicknessCell => layerThicknessCellWetDry
       ! need predicted transport velocity to limit drying flux
       !$omp parallel
-      !$omp do schedule(runtime) private(i, iEdge, k, divOutFlux, layerThickness)
+      !$omp do schedule(runtime) private(i, iEdge, k, divOutFlux)
       do iCell = 1, nCellsAll
         do k = minLevelCell(iCell), maxLevelCell(iCell)
           divOutFlux = 0.0_RKIND
-          layerThickness = min(layerThicknessProvis(k, iCell), layerThicknessCur(k, iCell))
           do i = 1, nEdgesOnCell(iCell)
             iEdge = edgesOnCell(i, iCell)
             if (k <= maxLevelEdgeTop(iEdge) .and. k >= minLevelEdgeBot(iEdge)) then
@@ -382,15 +407,148 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness
               end if
             end if
           end do
-          layerThickness = layerThickness + dt * divOutFlux
-
-          call ocn_wetting_velocity_factor_on_cell_edges(wettingVelocityFactor, layerThickness, normalVelocity, iCell, k)
 
+          layerThicknessCell(k, iCell) = min(layerThicknessProvis(k, iCell), layerThicknessCur(k, iCell)) + dt * divOutFlux
+          if ( .not. config_use_ssh_gradient_wetting_drying ) then
+             call ocn_wetting_velocity_factor_on_cell_edges( &
+                      wettingVelocityFactor, layerThicknessCell(k, iCell), normalVelocity, iCell, k)
+          end if
         end do
       end do
       !$omp end do
       !$omp end parallel
 
+      call mpas_dmpar_field_halo_exch(domain, 'layerThicknessCellWetDry')
+
+
+      ! This wetting and drying implementation is based on the one above.
+      ! However, it has special considerations for edges with normal velocities
+      ! that have been previously zeroed to limit drying (or are zero initially),
+      ! but have ssh gradients that would cause flow into the current cell,
+      ! i.e., the neighboring cell's ssh is higher.
+      ! There are two cases that require special treatment:
+      !  1) In the case that the neighbor cell is predicted to be wet, flow should
+      !     be allowed into the current cell.
+      !  2) In the case that the neighor cell is predicted to be dry, flow should
+      !     be prevented completely (and not ramped).
+      ! Currently, this only works for single layer configurations with the 'ssh_gradient'
+      ! pressure gradient option
+      if (config_use_ssh_gradient_wetting_drying ) then
+        ! Compute predicted ssh
+        sshCell => sshCellWetDry
+        sshCell = 0.0_RKIND
+
+        wettingVelocityFactor = 0.0_RKIND
+
+        if ( config_use_subgrid_wetting_drying ) then
+          do iCell = 1, nCellsAll
+            k = 1
+            call ocn_subgrid_ssh_lookup(layerThicknessCell(k, iCell), &
+                                        subgridWetVolumeCellTable(:,iCell), &
+                                        subgridSshCellTableRange(:,iCell), &
+                                        bottomDepth(iCell), &
+                                        subgridCellBathymetryMin(iCell), &
+                                        sshCell(iCell))
+          enddo
+        else
+          !$omp parallel
+          !$omp do schedule(runtime) private(k, columnThickness)
+          do iCell = 1, nCellsAll
+            columnThickness = 0.0_RKIND
+            do k = minLevelCell(iCell), maxLevelCell(iCell)
+              columnThickness = columnThickness + layerThicknessCell(k, iCell)
+            enddo
+            sshCell(iCell) = columnThickness - bottomDepth(iCell)
+          enddo
+          !$omp end do
+          !$omp end parallel
+        end if
+
+        !$omp parallel
+        !$omp do schedule(runtime) private(k, i, iEdge, cellDummy, cellNei, cellCur, sshCur, sshNei)
+        do iCell = 1, nCellsAll
+          do i = 1, nEdgesOnCell(iCell)
+            iEdge = edgesOnCell(i, iCell)
+            do k = minLevelCell(iCell), maxLevelCell(iCell)
+              if (k <= maxLevelEdgeTop(iEdge) .and. k >= minLevelEdgeBot(iEdge) .and. &
+                 layerThicknessCell(k, iCell) <= hCrit ) then
+                ! If edge velocity has been zeroed, check neighbor and currnt cell ssh.
+                ! Allow flow into current cell if neighbor cell is wet and the ssh gradient
+                ! would cause flow from neighbor cell into current cell.
+                if ( normalVelocity(k,iEdge)*edgeSignOnCell(i, iCell) < 0.0_RKIND ) then
+
+                  wettingVelocityFactor(k,iEdge) = 1.0_RKIND
+
+                elseif ( normalVelocity(k,iEdge)*edgeSignOnCell(i, iCell) == 0.0_RKIND ) then
+
+                  cellDummy(1:2) = cellsOnEdge(1:2,iEdge)
+                  cellNei = merge( cellDummy(2), cellDummy(1), iCell == cellDummy(1) )
+                  cellCur = iCell
+                  if ( cellNei > nCellsAll ) cellNei = cellCur
+                  sshCur = sshCell(cellCur)
+                  sshNei = sshCell(cellNei)
+
+                  wettingVelocityFactor(k,iEdge) = 1.0_RKIND
+
+                  ! if the neigbor cell is anticipated wet
+                  if ( layerThicknessCell(k, cellNei) >  hCrit ) then
+                    if ( sshCur < sshNei ) then
+                       wettingVelocityFactor(k,iEdge) = 0.0_RKIND
+                    end if
+                  end if
+
+                end if ! velocity check
+              endif ! k check and thickness check
+            enddo ! k
+          enddo ! i
+        end do ! iCell
+        !$omp end do
+        !$omp end parallel
+
+        !$omp parallel
+        !$omp do schedule(runtime) private(k, i, iEdge, cellDummy, cellNei, cellCur, sshCur, sshNei)
+        do iCell = 1, nCellsAll
+          do i = 1, nEdgesOnCell(iCell)
+            iEdge = edgesOnCell(i, iCell)
+            do k = minLevelCell(iCell), maxLevelCell(iCell)
+              if (k <= maxLevelEdgeTop(iEdge) .and. k >= minLevelEdgeBot(iEdge) .and. &
+                  config_zero_drying_velocity_ramp .and. &
+                  (layerThicknessCell(k, iCell) > hCrit) .and. &
+                  (layerThicknessCell(k, iCell) <= hRampMax) ) then
+
+                if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then
+                  wettingVelocityFactor(k, iEdge) = 1.0_RKIND - &
+                      tanh(50.0_RKIND * (layerThicknessCell(k, iCell) - hRampMin)/hRampMin)
+                endif
+
+                ! If edge velocity has been zeroed, check neighbor and current cell ssh.
+                ! Prevent flow out of cell if neighbor cell is dry and the ssh gradient
+                ! would cause flow from neighbor cell into current cell.
+                if ( normalVelocity(k, iEdge) * edgeSignOnCell(i,iCell) == 0.0_RKIND ) then
+
+                  cellDummy(1:2) = cellsOnEdge(1:2,iEdge)
+                  cellNei = merge( cellDummy(2), cellDummy(1), iCell == cellDummy(1) )
+                  cellCur = iCell
+                  if ( cellNei > nCellsAll ) cellNei = cellCur
+                  sshCur = sshCell(cellCur)
+                  sshNei = sshCell(cellNei)
+
+                  ! if the neigbor cell is anticipated dry
+                  if ( layerThicknessCell(k, cellNei) <=  hCrit ) then
+                     if ( sshCur < sshNei ) then
+                         wettingVelocityFactor(k,iEdge) = 1.0_RKIND
+                     end if
+                  end if
+                end if ! velocity check
+              end if ! thickness check
+            enddo ! k
+          enddo ! iEdge
+        end do ! iCell
+        !$omp end do
+        !$omp end parallel
+
+      end if
+
    end subroutine ocn_wetting_drying_wettingVelocity !}}}
 
 !***********************************************************************