Add subgrid scale correction scheme and init mode test cases

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

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

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 #
 ################################################

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>

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"

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">
 		<!-- **********************************************************************

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
       !

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)

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

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)
 

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"