Merge branch 'cbegeman/ocn/add-wetting-drying-ramp-feature' into next…

… (PR #5590)

Add ocean wetting-and-drying ramp feature

This PR enhances the existing wetting-and-drying algorithm using an
approach from O'Dea et al. 2020. Instead of zeroing out normalVelocity
and normalVelocity tendencies in cells where the projected
layerThickness drops below the user-defined minimum, this method ramps
down the normalVelocity and its tendencies between a range of
layerThicknesses.

Wetting-and-drying is currently only used in standalone ocean runs.

[BFB]

components/mpas-ocean/src/Registry.xml

@@ -1121,14 +1121,26 @@
 					possible_values=".true. or .false."
 		/>
 		<nml_option name="config_drying_min_cell_height" type="real" default_value="1.0e-3" units="m"
-					description="Minimum allowable cell height under drying.  Cell to be kept wet to at least this thickness."
+					description="Minimum allowable cell height under drying.  Cell to be kept wet to at least this thickness. When ramp is applied this is the min edge height"
 					possible_values="any positive real, typically 1.0e-3"
 		/>
 		<nml_option name="config_zero_drying_velocity" type="logical" default_value=".false."
 					description="If true, just zero out velocity that is contributing to drying for cell that is drying.  This option can be used to estimate acceptable minimum thicknesses for a run."
 					possible_values=".true. or .false."
 		/>
-		<nml_option name="config_verify_not_dry" type="logical" default_value=".false."
+		<nml_option name="config_zero_drying_velocity_ramp" type="logical" default_value=".false."
+					description="If true, ramp velocities and tendencies to zero rather than applying a simple on/off switch."
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_zero_drying_velocity_ramp_hmin" type="real" default_value="1e-3"
+					description="Minimum layer thickness at which velocities and tendencies are ramped toward zero. Recommended value equal to config_drying_min_cell_height."
+					possible_values="Any positive real"
+		/>
+		<nml_option name="config_zero_drying_velocity_ramp_hmax" type="real" default_value="2e-3"
+					description="Maximum layer thickness at which velocities and tendencies are ramped toward zero. Recommended values between 2x and 10x config_drying_min_cell_height."
+					possible_values="Any positive real"
+		/>
+		<nml_option name="config_verify_not_dry" type="logical" default_value=".false." units="unitless"
 					description="If true, verify that cells are at least config_min_cell_height thick."
 					possible_values=".true. or .false."
 		/>

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration.F

@@ -218,7 +218,7 @@ subroutine ocn_timestep_init(domain, dt, err)!{{{
 
       case ('RK4')
          timeIntegratorChoice = timeIntRK4
-         ! No init routine needs to be run for RK4
+         call ocn_time_integration_rk4_init(domain)
 
       case default
 

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F

@@ -61,6 +61,7 @@ module ocn_time_integration_rk4
    !--------------------------------------------------------------------
 
    public :: ocn_time_integrator_rk4
+   public :: ocn_time_integration_rk4_init
 
    contains
 
@@ -669,7 +670,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
          end if
 
          if (config_disable_tr_all_tend) then
-            num_tracers = size(activeTracersNew, dim=2)
+            num_tracers = size(activeTracersNew, dim=1)
             do iTracer = 1, num_tracers
                !$omp parallel
                !$omp do schedule(runtime)
@@ -1702,6 +1703,51 @@ subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{
 
    end subroutine ocn_time_integrator_rk4_cleanup!}}}
 
+!***********************************************************************
+!
+!  routine ocn_time_integration_rk4_init
+!
+!> \brief   Initialize RK4 time stepping within ocean
+!> \author  Carolyn Begeman
+!> \date    April 2023
+!> \details
+!>  This routine checks config options for RK4 integrator
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_time_integration_rk4_init(domain)!{{{
+
+      !-----------------------------------------------------------------
+      ! input/output variables
+      !-----------------------------------------------------------------
+      type (domain_type), intent(in) :: &
+         domain      !< [inout] data structure containing most variables
+      !-----------------------------------------------------------------
+      ! local variables
+      !-----------------------------------------------------------------
+      type (block_type), pointer :: block
+      type (mpas_pool_type), pointer :: meshPool
+      logical, pointer :: config_use_debugTracers
+      integer, pointer :: nVertLevels
+      ! End preamble
+      !-----------------------------------------------------------------
+      ! Begin code
+
+      call mpas_pool_get_config(domain % configs, 'config_use_debugTracers', config_use_debugTracers)
+      block => domain % blocklist
+      do while (associated(block))
+         call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
+         call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
+         if (config_use_debugTracers .and. nVertLevels == 1) then
+            call mpas_log_write('Debug tracers cannot be used in a ' &
+               // 'single layer case. Consider setting ' &
+               // 'config_use_debugTracers to .false.', MPAS_LOG_CRIT)
+         endif
+         block => block % next
+      end do
+
+   end subroutine ocn_time_integration_rk4_init
+
 end module ocn_time_integration_rk4
 
 ! vim: foldmethod=marker

components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F

@@ -261,7 +261,8 @@ 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 sure tendency doesn't dry cells
+     ! 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, &
                                              normalTransportVelocity, rkSubstepWeight, wettingVelocityFactor, err)
@@ -274,8 +275,10 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying
          do k = minLevelEdgeTop(iEdge), maxLevelEdgeBot(iEdge)
 
            if (abs(wettingVelocityFactor(k, iEdge)) > 0.0_RKIND) then
-             normalTransportVelocity(k, iEdge) = 0.0_RKIND
-             normalVelocity(k, iEdge) = 0.0_RKIND
+             normalTransportVelocity(k, iEdge) = (1.0_RKIND - &
+               wettingVelocityFactor(k, iEdge)) * normalTransportVelocity(k, iEdge)
+             normalVelocity(k, iEdge) = (1.0_RKIND - &
+               wettingVelocityFactor(k, iEdge)) * normalVelocity(k, iEdge)
            end if
 
          end do
@@ -347,20 +350,24 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness
       !
       !-----------------------------------------------------------------
 
-      integer :: iEdge, iCell, k, i
+      integer :: cell1, cell2, iEdge, iCell, k, i
 
-      real (kind=RKIND) :: divOutFlux
+      real (kind=RKIND) :: divFlux, divOutFlux
       real (kind=RKIND) :: layerThickness
+      real (kind=RKIND) :: hCrit, hRampMin, hEdgeTotal
+
+      character (len=100) :: log_string
 
       err = 0
 
-      if (.not. config_zero_drying_velocity ) return
+      hCrit = config_drying_min_cell_height
+
+      if (.not. config_zero_drying_velocity) return
 
       ! need predicted transport velocity to limit drying flux
       !$omp parallel
       !$omp do schedule(runtime) private(i, iEdge, k, divOutFlux, layerThickness)
       do iCell = 1, nCellsAll
-        ! can switch with maxLevelEdgeBot(iEdge)
         do k = minLevelCell(iCell), maxLevelCell(iCell)
           divOutFlux = 0.0_RKIND
           layerThickness = min(layerThicknessProvis(k, iCell), layerThicknessCur(k, iCell))
@@ -369,27 +376,44 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness
             if (k <= maxLevelEdgeTop(iEdge) .and. k >= minLevelEdgeBot(iEdge)) then
               ! only consider divergence flux leaving the cell
               if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) < 0.0_RKIND ) then
-                divOutFlux = divOutFlux &
-                             + normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) * &
-                               layerThickEdgeFlux(k, iEdge) * dvEdge(iEdge) * &
-                               invAreaCell(iCell)
+                divOutFlux = divOutFlux + &
+                             normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) * &
+                             layerThickEdgeFlux(k, iEdge) * dvEdge(iEdge) * &
+                             invAreaCell(iCell)
               end if
             end if
           end do
+          layerThickness = layerThickness + dt * divOutFlux
 
-          ! if layer thickness is too small, limit divergence flux outwards with opposite velocity
-          if ((layerThickness + dt * divOutFlux ) <= &
-              (config_drying_min_cell_height + config_drying_safety_height)) then
-
+          ! if layer thickness is too small, limit divergence flux outwards with
+          ! opposite velocity
+          if (layerThickness <= &
+              hCrit + config_drying_safety_height) then
             do i = 1, nEdgesOnCell(iCell)
               iEdge = edgesOnCell(i, iCell)
               if (k <= maxLevelEdgeBot(iEdge) .and. k >= minLevelEdgeTop(iEdge)) then
                 if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then
-                  ! just go with simple boolean approach for zero wetting velocity for debugging purposes
                   wettingVelocityFactor(k, iEdge) = 1.0_RKIND
                 end if
               end if
             end do
+          elseif (config_zero_drying_velocity_ramp .and. &
+                 (layerThickness > &
+                 hCrit + config_drying_safety_height) .and. &
+                 (layerThickness <= config_zero_drying_velocity_ramp_hmax)) then
+
+            hRampMin = config_zero_drying_velocity_ramp_hmin
+            ! Following O'Dea et al. (2021), if total upwinded wct is less than
+            ! 2*critical thickness, apply damping at each edge
+            do i = 1, nEdgesOnCell(iCell)
+              iEdge = edgesOnCell(i, iCell)
+              if (k <= maxLevelEdgeBot(iEdge) .and. k >= minLevelEdgeTop(iEdge)) then
+                if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then
+                  wettingVelocityFactor(k, iEdge) = 1.0_RKIND - &
+                    tanh(50.0_RKIND * (layerThickness - hRampMin)/hRampMin)
+                end if
+              end if
+            end do
 
           end if