Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Clean-up wetting and drying routines #14

Closed
wants to merge 7 commits into from
Closed

Clean-up wetting and drying routines #14

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
029e33e
Use ocn_mesh in ocn_wetting_drying
cbegeman Mar 1, 2022
ada49ff
Change wettingVelocity to wettingVelocityFactor
cbegeman Mar 1, 2022
b8b98fa
Clean up wetting and drying
cbegeman Mar 1, 2022
666371a
Move conditionals for efficiency
cbegeman Apr 18, 2022
2e6a16d
Fixup Registry for consistency with CF work
cbegeman Jul 13, 2022
f93a4b1
Aesthetic change to ocn_wetting_drying
cbegeman Jul 13, 2022
8ff26b8
Bugfix: use dt not rk4 substep in wettingVelocity routine
cbegeman Sep 26, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions components/mpas-ocean/src/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2693,8 +2693,8 @@
description="horizontal velocity used to transport mass and tracers"
packages="forwardMode;analysisMode"
/>
<var name="wettingVelocity" type="real" dimensions="nVertLevels nEdges Time" units="m s^-1"
description="Velocity to prevent drying of cell."
<var name="wettingVelocityFactor" type="real" dimensions="nVertLevels nEdges Time"
description="Scaling factor for normalVelocity and its tendency to prevent drying of cell between 0 (no scaling) and 1 (normalVelocity and tendency set to 0)."
packages="forwardMode"
/>
<var name="vertAleTransportTop" type="real" dimensions="nVertLevelsP1 nCells Time" units="m s^-1"
Expand Down
6 changes: 3 additions & 3 deletions components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1231,7 +1231,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight,
do k = 1, nVertLevels
normalVelocityProvis(k, iEdge) = edgeMask(k, iEdge)*(normalVelocityCur(k, iEdge) + rkSubstepWeight &
* normalVelocityTend(k, iEdge))
normalVelocityProvis(k, iEdge) = normalVelocityProvis(k, iEdge) * (1.0_RKIND - wettingVelocity(k, iEdge))
normalVelocityProvis(k, iEdge) = normalVelocityProvis(k, iEdge) * (1.0_RKIND - wettingVelocityFactor(k, iEdge))
end do
end do
!$omp end do
Expand Down Expand Up @@ -1474,12 +1474,12 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{

! the loop below is on all vertlevels but we used edgemask to set
! velocities to zero below topo (esp. edges between topo and ocean)
!$omp do schedule(runtime)
!$omp do schedule(runtime) private(k)
do iEdge = 1, nEdges
do k = 1, nVertLevels
normalVelocityNew(k, iEdge) = edgeMask(k,iEdge) &
*(normalVelocityNew(k, iEdge) + rkWeight * normalVelocityTend(k, iEdge))
normalVelocityNew(k, iEdge) = normalVelocityNew(k, iEdge) * (1.0_RKIND - wettingVelocity(k, iEdge))
normalVelocityNew(k, iEdge) = normalVelocityNew(k, iEdge) * (1.0_RKIND - wettingVelocityFactor(k, iEdge))
end do
end do
!$omp end do
Expand Down
2 changes: 1 addition & 1 deletion components/mpas-ocean/src/shared/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -216,7 +216,7 @@ mpas_ocn_framework_forcing.o:

mpas_ocn_time_varying_forcing.o: mpas_ocn_framework_forcing.o mpas_ocn_diagnostics_variables.o mpas_ocn_constants.o mpas_ocn_config.o

mpas_ocn_wetting_drying.o: mpas_ocn_diagnostics.o mpas_ocn_gm.o mpas_ocn_diagnostics_variables.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_gm.o
mpas_ocn_wetting_drying.o: mpas_ocn_diagnostics.o mpas_ocn_gm.o mpas_ocn_diagnostics_variables.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_gm.o mpas_ocn_mesh.o

mpas_ocn_tidal_potential_forcing.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o

Expand Down
16 changes: 8 additions & 8 deletions components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ module ocn_diagnostics_variables
real (kind=RKIND), dimension(:,:), pointer :: viscosity
real (kind=RKIND), dimension(:,:), pointer :: vertViscTopOfEdge
real (kind=RKIND), dimension(:,:), pointer :: vertViscTopOfCell
real (kind=RKIND), dimension(:,:), pointer :: wettingVelocity
real (kind=RKIND), dimension(:,:), pointer :: wettingVelocityFactor
type (field2DReal), pointer :: wettingVelocityField

real (kind=RKIND), dimension(:,:), pointer :: vertDiffTopOfCell
Expand Down Expand Up @@ -242,9 +242,9 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
end if

if ( trim(config_ocean_run_mode) == 'forward' ) then
call mpas_pool_get_array(diagnosticsPool, 'wettingVelocity', &
wettingVelocity)
call mpas_pool_get_field(diagnosticsPool, 'wettingVelocity', &
call mpas_pool_get_array(diagnosticsPool, 'wettingVelocityFactor', &
wettingVelocityFactor)
call mpas_pool_get_field(diagnosticsPool, 'wettingVelocityFactor', &
wettingVelocityField)
end if

Expand Down Expand Up @@ -763,8 +763,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
!$acc )
end if
if ( trim(config_ocean_run_mode) == 'forward' ) then
!$acc enter data create( &
!$acc wettingVelocity &
!$acc enter data copyin( &
!$acc wettingVelocityFactor &
!$acc )
end if
if (config_use_GM.or.config_use_Redi) then
Expand Down Expand Up @@ -1012,7 +1012,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
end if
if ( trim(config_ocean_run_mode) == 'forward' ) then
!$acc exit data delete( &
!$acc wettingVelocity &
!$acc wettingVelocityFactor &
!$acc )
end if
if (config_use_GM.or.config_use_Redi) then
Expand Down Expand Up @@ -1216,7 +1216,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
normalizedRelativeVorticityCell)
end if
if ( trim(config_ocean_run_mode) == 'forward' ) then
nullify(wettingVelocity)
nullify(wettingVelocityFactor)
end if
if (config_use_GM.or.config_use_Redi) then
nullify(RediKappa, &
Expand Down
7 changes: 3 additions & 4 deletions components/mpas-ocean/src/shared/mpas_ocn_tendency.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -351,7 +351,7 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, &
!$acc density, normRelVortEdge, montgomeryPotential, pressure, &
!$acc thermExpCoeff, salineContractCoeff, tangentialVelocity, &
!$acc layerThickEdgeFlux, kineticEnergyCell, sfcFlxAttCoeff, potentialDensity, &
!$acc vertAleTransportTop, vertViscTopOfEdge, wettingVelocity)
!$acc vertAleTransportTop, vertViscTopOfEdge, wettingVelocityFactor)
!$acc enter data &
!$acc copyin(tendVel, sfcStress, sfcStressMag, &
!$acc ssh, normalVelocity, &
Expand Down Expand Up @@ -452,15 +452,14 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, &

#ifdef MPAS_OPENACC
!$acc parallel loop collapse(2) &
!$acc present(tendVel, wettingVelocity)
!$acc present(tendVel, wettingVelocityFactor)
#else
!$omp parallel
!$omp do schedule(runtime) private(k)
#endif
do iEdge = 1, nEdgesAll
do k=1,nVertLevels
tendVel(k,iEdge) = tendVel(k,iEdge)* &
(1.0_RKIND - wettingVelocity(k,iEdge))
tendVel(k,iEdge) = tendVel(k,iEdge) * (1.0_RKIND - wettingVelocityFactor(k,iEdge))
end do
end do
#ifndef MPAS_OPENACC
Expand Down
125 changes: 39 additions & 86 deletions components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ module ocn_wetting_drying
use ocn_diagnostics
use ocn_diagnostics_variables
use ocn_gm
use ocn_mesh

implicit none
private
Expand Down Expand Up @@ -108,16 +109,13 @@ subroutine ocn_wetting_drying_verify( block , minHeight, err)!{{{
!-----------------------------------------------------------------

type (mpas_pool_type), pointer :: statePool, meshPool, tendPool
integer, dimension(:), pointer :: minLevelCell, maxLevelCell
real (kind=RKIND), dimension(:), pointer :: sshSubcycleNew
real (kind=RKIND), dimension(:), pointer :: bottomDepth
integer, pointer :: nCellsSolve
integer :: iCell, k
integer :: debugUnit
real (kind=RKIND), dimension(:,:), pointer :: layerThicknessCur
real (kind=RKIND), dimension(:,:), pointer :: layerThicknessNew
real (kind=RKIND), dimension(:,:), pointer :: layerThicknessTend
real (kind=RKIND), dimension(:), pointer :: lonCell, latCell
real (kind=RKIND) :: minThickness, layerThick
character (len=StrKIND) :: debugFilename

Expand All @@ -134,13 +132,8 @@ subroutine ocn_wetting_drying_verify( block , minHeight, err)!{{{
call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, timeLevel=1)
call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessNew, timeLevel=2)
call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve)
call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell)
call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
call mpas_pool_get_array(tendPool, 'layerThickness', layerThicknessTend)
call mpas_pool_get_array(statePool, 'sshSubcycle', sshSubcycleNew, 2)
call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth)
call mpas_pool_get_array(meshPool, 'lonCell', lonCell)
call mpas_pool_get_array(meshPool, 'latCell', latCell)

err = 0

Expand All @@ -151,8 +144,11 @@ subroutine ocn_wetting_drying_verify( block , minHeight, err)!{{{
do iCell = 1, nCellsSolve
do k = minLevelCell(iCell), maxLevelCell(iCell)
! use ssh as a proxy too for baroclinic mode
if (trim(config_time_integrator) == 'split_explicit' .or. trim(config_time_integrator) == 'split_implicit') then
layerThick = min(layerThicknessNew(k, iCell), (sshSubcycleNew(iCell)+bottomDepth(iCell))/maxLevelCell(iCell))
! Note: wetting-drying currently not supported for either of these time integration methods
if (trim(config_time_integrator) == 'split_explicit' .or. &
trim(config_time_integrator) == 'split_implicit') then
layerThick = min(layerThicknessNew(k, iCell), &
(sshSubcycleNew(iCell) + bottomDepth(iCell))/maxLevelCell(iCell))
else
layerThick = layerThicknessNew(k, iCell)
end if
Expand Down Expand Up @@ -237,23 +233,17 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying


type (mpas_pool_type), pointer :: tendPool
type (mpas_pool_type), pointer :: meshPool
type (mpas_pool_type), pointer :: statePool
type (mpas_pool_type), pointer :: provisStatePool
real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur
real (kind=RKIND), dimension(:, :), pointer :: layerThicknessProvis
real (kind=RKIND), dimension(:, :), pointer :: normalVelocity

integer, dimension(:), pointer :: minLevelEdgeTop
integer, dimension(:), pointer :: maxLevelEdgeTop
integer, dimension(:), pointer :: maxLevelEdgeBot
integer, pointer :: nEdges
integer :: iEdge, k

err = 0

call mpas_pool_get_subpool(block % structs, 'tend', tendPool)
call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
call mpas_pool_get_subpool(block % structs, 'state', statePool)
call mpas_pool_get_subpool(block % structs, 'provis_state', provisStatePool)

Expand All @@ -262,47 +252,37 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying
call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1)
call mpas_pool_get_array(provisStatePool, 'layerThickness', layerThicknessProvis, 1)

call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges)
call mpas_pool_get_array(meshPool, 'minLevelEdgeTop', minLevelEdgeTop)
call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop)
call mpas_pool_get_array(meshPool, 'maxLevelEdgeBot', maxLevelEdgeBot)

!$omp parallel
!$omp do schedule(runtime)
do iEdge = 1, nEdges
wettingVelocity(:, iEdge) = 0.0_RKIND
do iEdge = 1, nEdgesAll
wettingVelocityFactor(:, iEdge) = 0.0_RKIND
end do
!$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

call ocn_wetting_drying_wettingVelocity(meshPool, layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
normalTransportVelocity, rkSubstepWeight, wettingVelocity, err)
call ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
normalTransportVelocity, dt, wettingVelocityFactor, err)

! prevent drying from happening with selective wettingVelocity
!$omp parallel
!$omp do schedule(runtime) private(k)
do iEdge = 1, nEdges
do k = minLevelEdgeTop(iEdge), maxLevelEdgeBot(iEdge)
if (abs(normalTransportVelocity(k,iEdge) + wettingVelocity(k,iEdge)) < eps) then
! prevent spurious flux for close to zero values
normalTransportVelocity(k, iEdge) = 0.0_RKIND
normalVelocity(k, iEdge) = 0.0_RKIND
else if (abs(normalTransportVelocity(k,iEdge) + wettingVelocity(k,iEdge)) <= abs(normalTransportVelocity(k,iEdge))) then
normalTransportVelocity(k, iEdge) = normalTransportVelocity(k, iEdge) + wettingVelocity(k, iEdge)
normalVelocity(k, iEdge) = normalVelocity(k, iEdge) + wettingVelocity(k, iEdge)
end if

if (abs(wettingVelocity(k, iEdge)) > 0.0_RKIND .and. config_zero_drying_velocity) then
normalTransportVelocity(k, iEdge) = 0.0_RKIND
normalVelocity(k, iEdge) = 0.0_RKIND
end if
! prevent drying from happening with selective wettingVelocityFactor
if (config_zero_drying_velocity) then
!$omp parallel
!$omp do schedule(runtime) private(k)
do iEdge = 1, nEdgesAll
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
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_prevent_drying_rk4 !}}}

Expand All @@ -319,19 +299,16 @@ end subroutine ocn_prevent_drying_rk4 !}}}
!> to prevent cells from drying.
!
!-----------------------------------------------------------------------
subroutine ocn_wetting_drying_wettingVelocity(meshPool, layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &

normalVelocity, dt, wettingVelocity, err)!{{{
normalVelocity, dt, wettingVelocityFactor, err)!{{{

!-----------------------------------------------------------------
!
! input variables
!
!-----------------------------------------------------------------

type (mpas_pool_type), intent(in) :: &
meshPool !< Input: horizonal mesh information

real (kind=RKIND), dimension(:,:), intent(in) :: &
layerThicknessCur !< Input: layer thickness at old time

Expand All @@ -354,7 +331,7 @@ subroutine ocn_wetting_drying_wettingVelocity(meshPool, layerThickEdgeFlux, laye
!-----------------------------------------------------------------

real (kind=RKIND), dimension(:,:), intent(inout) :: &
wettingVelocity !< Input/Output: velocity wettingVelocityency
wettingVelocityFactor !< Input/Output: velocity wettingVelocityFactor

!-----------------------------------------------------------------
!
Expand All @@ -371,38 +348,18 @@ subroutine ocn_wetting_drying_wettingVelocity(meshPool, layerThickEdgeFlux, laye
!-----------------------------------------------------------------

integer :: iEdge, iCell, k, i
integer, pointer :: nVertLevels, nCells
integer, dimension(:), pointer :: nEdgesOnCell
integer, dimension(:), pointer :: minLevelCell, maxLevelCell
integer, dimension(:), pointer :: minLevelEdgeTop, maxLevelEdgeBot
integer, dimension(:,:), pointer :: edgesOnCell
integer, dimension(:,:), pointer :: edgeSignOnCell

real (kind=RKIND) :: divOutFlux
real (kind=RKIND) :: invAreaCell
real (kind=RKIND) :: layerThickness
real (kind=RKIND), dimension(:), pointer :: dvEdge
real (kind=RKIND), dimension(:), pointer :: areaCell

err = 0

call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell)
call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell)
call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell)
call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
call mpas_pool_get_array(meshPool, 'minLevelEdgeTop', minLevelEdgeTop)
call mpas_pool_get_array(meshPool, 'maxLevelEdgeBot', maxLevelEdgeBot)
call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge)
call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
if (.not. config_zero_drying_velocity ) return

! need predicted transport velocity to limit drying flux
!$omp parallel
!$omp do schedule(runtime) private(invAreaCell, i, iEdge, k, divOutFlux, layerThickness)
do iCell = 1, nCells
invAreaCell = 1.0_RKIND / areaCell(iCell)
!$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
Expand All @@ -412,32 +369,28 @@ subroutine ocn_wetting_drying_wettingVelocity(meshPool, layerThickEdgeFlux, laye
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
divOutFlux = divOutFlux &
+ normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) * &
layerThickEdgeFlux(k, iEdge) * dvEdge(iEdge) * &
invAreaCell(iCell)
end if
end if
end do

! 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
! limit divOutFlux out of cell to keep it wet
divOutFlux = abs(divOutFlux)
divOutFlux = (layerThickness - (config_drying_min_cell_height + eps)) / (dt*divOutFlux + eps)
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
! each outgoing velocity is penalized (but not the incoming, wetting velocities)
! square the fractional term to make values near zero go to zero much quicker (to prevent threshold from being hit)
wettingVelocity(k, iEdge) = - (min(max(0.0_RKIND, 1.0_RKIND - (divOutFlux*divOutFlux)), 1.0_RKIND)) * normalVelocity(k, iEdge)
! just go with simple boolean approach for zero wetting velocity for debugging purposes
if (config_zero_drying_velocity) then
wettingVelocity(k, iEdge) = 1.0_RKIND
end if
wettingVelocityFactor(k, iEdge) = 1.0_RKIND
end if
end if
end do

end if

end do
Expand Down