diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index fdf915646623..0d8d366e09b5 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1089,6 +1089,10 @@ 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." /> + h_min call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1) @@ -261,10 +263,12 @@ 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) + call ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, & + layerThicknessProvis, normalTransportVelocity, & + ssh, rkSubstepWeight, wettingVelocityFactor, err) ! prevent drying from happening with selective wettingVelocityFactor if (config_zero_drying_velocity) then @@ -274,8 +278,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 @@ -301,7 +307,7 @@ end subroutine ocn_prevent_drying_rk4 !}}} !----------------------------------------------------------------------- subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, & - normalVelocity, dt, wettingVelocityFactor, err)!{{{ + normalVelocity, ssh, dt, wettingVelocityFactor, err)!{{{ !----------------------------------------------------------------- ! @@ -321,6 +327,8 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness real (kind=RKIND), dimension(:,:), intent(in) :: & normalVelocity !< Input: transport + real (kind=RKIND), dimension(:), intent(in) :: ssh + real (kind=RKIND), intent(in) :: & dt !< Input: time step @@ -347,55 +355,86 @@ 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) :: layerThickness + real (kind=RKIND) :: hCrit, hEdgeTotal + + character (len=100) :: log_string err = 0 - if (.not. config_zero_drying_velocity ) return + 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)) - do i = 1, nEdgesOnCell(iCell) - iEdge = edgesOnCell(i, iCell) - if (k <= maxLevelEdgeBot(iEdge) .or. k >= minLevelEdgeTop(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) - end if - end if - end do + if (config_zero_drying_velocity_ramp) then + + ! Following O'Dea et al. (2021), if total upwinded wct is less than + ! 2*critical thickness, apply damping throughout the water column + + !$omp parallel + !$omp do schedule(runtime) private(cell1, cell2, iEdge, k, hEdgeTotal) + do iEdge = 1, nEdgesAll + + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + hEdgeTotal = 0.0_RKIND + hCrit = 0.0_RKIND + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + hEdgeTotal = hEdgeTotal + layerThickEdgeFlux(k, iEdge) + hCrit = hCrit + config_drying_min_cell_height + enddo - ! 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 (hEdgeTotal <= hCrit) then + wettingVelocityFactor(:, iEdge) = 1.0_RKIND + elseif (hEdgeTotal > hCrit .and. hEdgeTotal <= 2.0_RKIND*hCrit) then + wettingVelocityFactor(:, iEdge) = 1.0_RKIND - & + TANH(50.0_RKIND * (hEdgeTotal - hCrit)/hCrit) + endif + enddo + !$omp end do + !$omp end parallel + else + + ! need predicted transport velocity to limit drying flux + !$omp parallel + !$omp do schedule(runtime) private(i, iEdge, k, divOutFlux, layerThickness) + 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 <= maxLevelEdgeBot(iEdge) .or. 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 + ! 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) end if end if end do - end if + ! 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 + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + if (k <= maxLevelEdgeBot(iEdge) .or. k >= minLevelEdgeTop(iEdge)) then + if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then + wettingVelocityFactor(k, iEdge) = 1.0_RKIND + end if + end if + end do + end if + + end do end do - end do - !$omp end do - !$omp end parallel + !$omp end do + !$omp end parallel + end if end subroutine ocn_wetting_drying_wettingVelocity !}}}