From 241989bdeef97952724a65efe5304136169d74fc Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 18 Apr 2022 13:17:33 -0500 Subject: [PATCH 1/6] Add wettingVelocityFactor ramp as in O'Dea --- components/mpas-ocean/src/Registry.xml | 12 +++- .../src/shared/mpas_ocn_wetting_drying.F | 63 +++++++++++++------ 2 files changed, 55 insertions(+), 20 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 73e417b9303b..894107a37c20 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1129,14 +1129,22 @@ possible_values=".true. or .false." /> - + + 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..ed176258e9ef 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F @@ -235,6 +235,7 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying type (mpas_pool_type), pointer :: tendPool type (mpas_pool_type), pointer :: statePool type (mpas_pool_type), pointer :: provisStatePool + real (kind=RKIND), dimension(:), pointer :: ssh real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur real (kind=RKIND), dimension(:, :), pointer :: layerThicknessProvis real (kind=RKIND), dimension(:, :), pointer :: normalVelocity @@ -247,6 +248,7 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying call mpas_pool_get_subpool(block % structs, 'state', statePool) call mpas_pool_get_subpool(block % structs, 'provis_state', provisStatePool) + call mpas_pool_get_array(statePool, 'ssh', ssh, 1) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 1) ! use thickness at n because constraint is h_n + dt*T_h > h_min call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1) @@ -261,10 +263,11 @@ 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) + normalTransportVelocity, ssh, rkSubstepWeight, wettingVelocityFactor, err) ! prevent drying from happening with selective wettingVelocityFactor if (config_zero_drying_velocity) then @@ -274,8 +277,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 +306,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 +326,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,20 +354,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, 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 +380,43 @@ 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 + + ! 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) .or. 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 - hCrit)/hCrit) + end if + end if + end do end if From 9933271172b48f59fe7affc0b805f76e28552bf8 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 15 Sep 2022 16:42:56 -0500 Subject: [PATCH 2/6] Add minimum limit to ramp (different from hCrit) --- components/mpas-ocean/src/Registry.xml | 8 ++++++-- .../mpas-ocean/src/shared/mpas_ocn_wetting_drying.F | 11 +++++------ 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 894107a37c20..c95693a04bb0 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1140,9 +1140,13 @@ description="If true, ramp velocities and tendencies to zero rather than applying a simple on/off switch." possible_values=".true. or .false." /> + = 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 - hCrit)/hCrit) + tanh(50.0_RKIND * (layerThickness - hRampMin)/hRampMin) end if end if end do From 38968b3b61c3aa5a0a1527a50ea3278b70710787 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 9 Sep 2022 14:47:42 -0500 Subject: [PATCH 3/6] Change k-loop limits for wettingVelocityFactor --- components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 c58e2754c2b6..a288a6eaf829 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F @@ -409,7 +409,7 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness ! 2*critical thickness, apply damping at each edge do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) - if (k <= maxLevelEdgeBot(iEdge) .or. k >= minLevelEdgeTop(iEdge)) then + 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) From fc4375d7d19bc735a82973a9a497e74f0033c333 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 10 Apr 2023 16:59:42 -0500 Subject: [PATCH 4/6] Fix bug introduced by RK4 disable tracers PR#5519 --- .../mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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..b548bfc0f188 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 @@ -669,7 +669,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) From 770e7067d7db6f5d87d78b465d319d01b1ca923a Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 25 Apr 2023 09:17:39 -0700 Subject: [PATCH 5/6] fixup W&D ssh --- components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F | 2 -- 1 file changed, 2 deletions(-) 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 a288a6eaf829..aeb1dcdb4b7d 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F @@ -235,7 +235,6 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying type (mpas_pool_type), pointer :: tendPool type (mpas_pool_type), pointer :: statePool type (mpas_pool_type), pointer :: provisStatePool - real (kind=RKIND), dimension(:), pointer :: ssh real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur real (kind=RKIND), dimension(:, :), pointer :: layerThicknessProvis real (kind=RKIND), dimension(:, :), pointer :: normalVelocity @@ -248,7 +247,6 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying call mpas_pool_get_subpool(block % structs, 'state', statePool) call mpas_pool_get_subpool(block % structs, 'provis_state', provisStatePool) - call mpas_pool_get_array(statePool, 'ssh', ssh, 1) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 1) ! use thickness at n because constraint is h_n + dt*T_h > h_min call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1) From 73b90f7a392587840855cc6a05945a23e0209196 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 25 Apr 2023 19:42:13 -0700 Subject: [PATCH 6/6] Add error message for single_layer and debugTracers --- .../mode_forward/mpas_ocn_time_integration.F | 2 +- .../mpas_ocn_time_integration_rk4.F | 46 +++++++++++++++++++ 2 files changed, 47 insertions(+), 1 deletion(-) 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 b548bfc0f188..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 @@ -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