Buttermilk bay dev 2ndattempt

components/mpas-ocean/src/Registry.xml

@@ -1164,6 +1164,10 @@
 		<nml_option name="config_subgrid_table_levels_init" type="integer" default_value="10" units="unitless"
 					description="Number of ssh increments used in pre-computed subgrid lookup tables"
 					possible_values="Any positive integer"
+	        />
+		<nml_option name="config_use_modify_wetting_drying_implementation" type="logical" default_value=".false." units="unitless"
+					description="If true, use alternatove implementation of wetting and drying algorithm"
+					possible_values=".true. or .false."
 		/>
 	</nml_record>
 	<nml_record name="ocean_constants">
@@ -1550,7 +1554,6 @@
 		<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"/>
 		<package name="variableShortwave" description="This package includes variables required to compute spatially variable shortwave extinction coefficients"/>
-
 		<package name="splitTimeIntegrator" description="This package includes variables required for either the split or unsplit explicit time integrators."/>
 		<package name="semiImplicitTimePKG" description="This package includes variables required for split-implicit time integrators."/>
 		<package name="thicknessFilter" description="This package includes variables required for frequency filtered thickness."/>
@@ -2561,15 +2564,15 @@
 			 description="Pre-computed wet volume per unit area for each edge at given ssh value"
 			 packages="subgridWetDryPKG"
 		/>
-		<var name="subgridSshCellTableRange" type="real" dimensions="TWO nCells" units="m"
+		<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="TWO nEdges" units="m"
+		<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="TWO nVertices" units="m"
+		<var name="subgridSshVertexTableRange" type="real" dimensions="R3 nVertices" units="m"
 			 description="Range of ssh vaules used for vertex lookup table"
 			 packages="subgridWetDryPKG"
 		/>
@@ -2926,7 +2929,7 @@
 			 missing_value="FILLVAL" missing_value_mask="cellMask"
 			 packages="forwardMode;analysisMode"
 		/>
-		<var name="barotropicForcing" type="real" dimensions="nEdges Time" units="m s^-2"
+		<var name="barotropicForcing" type="real" dimensions="nEdges Time" uWCCM 2024 minisymposiumnits="m s^-2"
 			 description="Barotropic tendency computed from the baroclinic equations in stage 1 of the split-explicit algorithm."
 			 packages="forwardMode;analysisMode"
 		/>
@@ -3420,6 +3423,14 @@
 			description="gradient of sea surface height perturbation from self-attraction and loading"
 			packages="tidalPotentialForcingPKG"
 		/>
+		<var name="layerThicknessCellWetDry" type="real" dimensions="nCells Time" units="m"
+			description="intermediate estimate total layer thickness for wetting/drying purpose"
+			packages="forwardMode"
+		/>	
+		<var name="sshCellWetDry" type="real" dimensions="nCells Time" units="m"
+			description="intermediate estimate ssh for wetting/drying purpose"
+			packages="forwardMode"
+		/>
 	</var_struct>
 	<var_struct name="shortwave" time_levs="1">
 		<!-- **********************************************************************

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,10 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
 
       real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersCur, activeTracersNew
 
+      ! DW
+      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)
@@ -249,6 +256,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
 
@@ -480,10 +494,15 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
              ! compute wetting velocity to prevent drying of cell (sets up start of next iterate to not dry)
              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)
-                   block => block % next
-                 end do
+               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 = domain, &
+                        block = block, dt = dt, &
+                        rkSubstepWeight = rk_substep_weights(rk_step), &
+                        config_zero_drying_velocity = config_zero_drying_velocity, err = err)
+
+                  block => block % next
+               end do
                ! exchange fields for parallelization
                call mpas_pool_get_field(statePool, 'normalVelocity', normalVelocityField, 1)
                call mpas_dmpar_exch_halo_field(normalVelocityField)
@@ -708,6 +727,18 @@ 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
+
+               ! DW
+               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

@@ -31,7 +31,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

components/mpas-ocean/src/mode_init/Registry.xml

@@ -20,6 +20,8 @@
 #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"
 
 

components/mpas-ocean/src/mode_init/Registry_Buttermilk_Bay.xml

@@ -0,0 +1,75 @@
+	<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_Coriolis_parameter" type="real" default_value="0.0" units="1/s"
+			description="Coriolis paramter"
+			possible_values="Any real number"
+		/>
+	        <nml_option name="config_Buttermilk_bay_eta0" type="real" default_value="2.0" units="m"
+			description="Magntiude of elevation forcing boundary"
+			possible_values="Any real number"
+		/>
+		<nml_option name="config_Buttermilk_bay_omega" type="real" default_value="1.4544e-4" units="1/s"
+			description="Angular frequency of elevation forcing boundary"
+			possible_values="Any real number"
+		/>
+		<nml_option name="config_Buttermilk_bay_gravity" type="real" default_value="9.81" units="m/s^2"
+			description="Gravitational accerlation"
+			possible_values="Any real number"
+		/>
+	        <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_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_setup" type="logical" default_value="true" units="unitless"
+			description="logical flag controling the initial ssh and layerThichmess. standard setup is obtained by setting this namelist to .false."
+			possible_values=".true. or .false."
+		/>
+		<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_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>
+
+
+
+
+

components/mpas-ocean/src/mode_init/Registry_parabolic_bowl.xml

@@ -31,7 +31,11 @@
 		<nml_option name="config_parabolic_bowl_subgrid_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="Level of refinement to evaluate subgrid bathymety"
+			possible_values=".true. or .false."
+		/>
 	</nml_record>