diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index ff9d06c1c8bc..467d831ff299 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -1121,14 +1121,26 @@
possible_values=".true. or .false."
/>
-
+
+
+
diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration.F
index ba277921c60d..82ac93a7b129 100644
--- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration.F
+++ b/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
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 6aeb8274dea7..58755c871d12 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
@@ -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
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 75005d87c103..aeb1dcdb4b7d 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F
+++ b/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