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

Jump to bottom

Fixes MPAS-Ocean Richardson number calculation #5946

Merged
merged 7 commits into from
Oct 26, 2023
Merged

Fixes MPAS-Ocean Richardson number calculation #5946

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
4ce8596
Changes to richardson number calculation
vanroekel Sep 6, 2023
b853e01
Fixes halo on tangentialVelocity
vanroekel Sep 15, 2023
efacea4
Fixes to cvmix and domain decomp
vanroekel Sep 19, 2023
98540fe
Rolls back some extra halos
vanroekel Sep 19, 2023
9cd7f13
Adds halo exchange for RK4 and SI
vanroekel Sep 19, 2023
5c6f47c
add domain to RK4 subroutine
mark-petersen Sep 19, 2023
cfb1208
Add missing variable to omp directive
jonbob Oct 25, 2023
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
7 changes: 5 additions & 2 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 @@ -603,7 +603,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
! Rescale tracers
block => domain % blocklist
do while(associated(block))
call ocn_time_integrator_rk4_cleanup(block, dt, err)
call ocn_time_integrator_rk4_cleanup(domain, block, dt, err)

block => block % next
end do
Expand Down Expand Up @@ -1549,7 +1549,8 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{

end subroutine ocn_time_integrator_rk4_accumulate_update!}}}

subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{
subroutine ocn_time_integrator_rk4_cleanup(domain, block, dt, err)!{{{
type (domain_type), intent(inout) :: domain !< Input/Output: domain information
type (block_type), intent(in) :: block
real (kind=RKIND), intent(in) :: dt
integer, intent(out) :: err
Expand Down Expand Up @@ -1656,6 +1657,8 @@ subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{
endif
#endif
call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMeshPool, scratchPool, tracersPool, 2)
call mpas_dmpar_field_halo_exch(domain, 'tangentialVelocity')

#ifdef MPAS_OPENACC
!$acc update host(layerThickEdgeFlux, layerThickEdgeMean)
!$acc update host(relativeVorticity, circulation)
Expand Down
1 change: 1 addition & 0 deletions components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2840,6 +2840,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{
scratchPool, tracersPool, 2)

call mpas_dmpar_field_halo_exch(domain, 'surfaceFrictionVelocity')
call mpas_dmpar_field_halo_exch(domain, 'tangentialVelocity')

#ifdef MPAS_OPENACC
!$acc update host(layerThickEdgeFlux, layerThickEdgeMean)
Expand Down
2 changes: 1 addition & 1 deletion components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2402,7 +2402,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
scratchPool, tracersPool, 2)

call mpas_dmpar_field_halo_exch(domain, 'surfaceFrictionVelocity')

call mpas_dmpar_field_halo_exch(domain, 'tangentialVelocity')
#ifdef MPAS_OPENACC
!$acc update host(layerThickEdgeFlux, layerThickEdgeMean)
!$acc update host(relativeVorticity, circulation)
Expand Down
46 changes: 33 additions & 13 deletions components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1551,10 +1551,12 @@ subroutine ocn_diagnostic_solve_richardson(displacedDensity, density, zMid, &
!
!-----------------------------------------------------------------

real(kind=RKIND) :: coef, shearSquared, delU2, shearMean
integer :: iCell, k, i, iEdge
real(kind=RKIND) :: coef, shear, Ri_edge, dz_edge, shearSquared, delU2, shearMean
integer :: jCell, iCell, k, i, iEdge
integer :: nCells
real(kind=RKIND), dimension(:), allocatable :: RiTopOfCellNorm

allocate (RiTopOfCellNorm(nVertLevels))
!
! Brunt-Vaisala frequency (this has units of s^{-2})
!
Expand Down Expand Up @@ -1592,32 +1594,50 @@ subroutine ocn_diagnostic_solve_richardson(displacedDensity, density, zMid, &
!$acc parallel loop gang vector &
!$acc present(RiTopOfCell, maxLevelCell, nEdgesOnCell, edgesOnCell, normalVelocity, &
!$acc edgeAreaFractionOfCell, zMid, BruntVaisalaFreqTop, minLevelCell) &
!$acc private(k, i, iEdge, delU2, shearSquared, shearMean)
!$acc private(k, i, iEdge, delU2, shearMean, jCell, shear, dz_edge, Ri_edge, &
!$acc RiTopOfCellNorm)
#else
!$omp parallel
!$omp do schedule(runtime) private(k, shearSquared, i, iEdge, delU2, shearMean)
!$omp do schedule(runtime) private(k, i, iEdge, delU2, shearMean, jCell, shear, &
!$omp dz_edge, Ri_edge, RiTopOfCellNorm)
#endif
do iCell=1,nCells
RiTopOfCell(:,iCell) = 100.0_RKIND
do k = minLevelCell(iCell)+1, maxLevelCell(iCell)
shearSquared = 0.0_RKIND
RiTopOfCellNorm(:) = 1.0E-12_RKIND
do i = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(i, iCell)
delU2 = (normalVelocity(k-1,iEdge) - normalVelocity(k,iEdge))**2
shearSquared = shearSquared + edgeAreaFractionOfCell(i,iCell) * delU2
jCell = cellsOnCell(i,iCell)
do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge)
delU2 = (normalVelocity(k-1,iEdge) - normalVelocity(k,iEdge))**2 + &
( tangentialVelocity(k-1,iEdge) - tangentialVelocity(k,iEdge))**2

dz_edge = 0.5_RKIND * ( &
zMid(k-1,iCell) + zMid(k-1,jCell) - &
zMid(k-0,iCell) - zMid(k-0,jCell) )

shear = delU2 / dz_edge**2

Ri_edge = max(0.0_RKIND, 0.5_RKIND * ( &
BruntVaisalaFreqTop(k,iCell) + &
BruntVaisalaFreqTop(k,jCell) ) )/ (shear + 1.0E-12_RKIND)

RiTopOfCellNorm(k) = RiTopOfCellNorm(k) &
+ 0.25_RKIND * dcEdge(iEdge) * dvEdge(iEdge)

RiTopOfCell(k,iCell) = &
RiTopOfCell(k,iCell) + 0.25_RKIND * dcEdge(iEdge) * dvEdge(iEdge) * Ri_edge
enddo
! Note that the factor of two is from averaging dot product to cell center on a C-grid
shearMean = sqrt(2.0_RKIND*shearSquared )
shearMean = shearMean / (zMid(k-1,iCell) - zMid(k,iCell))
RiTopOfCell(k,iCell) = BruntVaisalaFreqTop(k,iCell) / (shearMean**2 + 1.0e-10_RKIND)
enddo
do k = minLevelCell(iCell), maxLevelCell(iCell)
RiTopOfCell(k,iCell) = RiTopOfCell(k,iCell) / RiTopOfCellNorm(k)
end do
RiTopOfCell(minLevelCell(iCell),iCell) = RiTopOfCell(minLevelCell(iCell)+1,iCell)
end do
#ifndef MPAS_OPENACC
!$omp end do
!$omp end parallel
#endif

deallocate(RiTopOfCellNorm)
end subroutine ocn_diagnostic_solve_richardson!}}}

!***********************************************************************
Expand Down
35 changes: 23 additions & 12 deletions components/mpas-ocean/src/shared/mpas_ocn_vmix_cvmix.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -159,12 +159,13 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim
integer, dimension(:), allocatable :: surfaceAverageIndex

real (kind=RKIND) :: x, bulkRichardsonNumberStop, sfc_layer_depth
real (kind=RKIND) :: normalVelocityAv, delU2, areaSum, blTemp
real (kind=RKIND) :: normalVelocityAv, delU2, areaSum, blTemp,tangentialVelocityAv
real (kind=RKIND) :: sigma
real (kind=RKIND), dimension(:), allocatable :: Nsqr_iface, turbulentScalarVelocityScale, &
deltaVelocitySquared, normalVelocitySum, &
potentialDensitySum, RiTemp, &
layerThicknessSum, layerThicknessEdgeSum
layerThicknessSum, layerThicknessEdgeSum, &
tangentialVelocitySum
real (kind=RKIND), dimension(:), allocatable, target :: RiSmoothed, BVFSmoothed, OBLDepths, interfaceForcings
logical :: bulkRichardsonFlag

Expand Down Expand Up @@ -308,6 +309,7 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim
allocate(RiTemp(nVertLevels+1))

allocate(normalVelocitySum(nVertLevels))
allocate(tangentialVelocitySum(nVertLevels))
allocate(potentialDensitySum(nVertLevels))
allocate(surfaceAverageIndex(nVertLevels))
allocate(deltaVelocitySquared(nVertLevels))
Expand All @@ -333,9 +335,10 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim
! !$omp private(k, RiSmoothed, RiTemp, nsmooth, BVFSmoothed, alphaAngle, &
! !$omp langmuirEnhancementFactor, bulkRichardsonNumberStop, bulkRichardsonFlag, &
! !$omp topIndex, kIndexOBL, deltaVelocitySquared, sigma, OBLDepths, interfaceForcings, &
! !$omp surfaceAverageIndex, sfc_layer_depth, kav, i, iEdge, &
! !$omp surfaceAverageIndex, sfc_layer_depth, kav, i, iEdge, tangentialVelocitySum, &
! !$omp normalVelocitySum, layerThicknessEdgeSum, normalVelocityAv, delU2, potentialDensitySum, &
! !$omp layerThicknessSum, blTemp, zonalWavenumberCoeff, meridionalWavenumberCoeff) &
! !$omp layerThicknessSum, blTemp, zonalWavenumberCoeff, meridionalWavenumberCoeff, &
! !$omp tangentialVelocityAv) &
! !$omp firstprivate(cvmix_variables, Nsqr_iface, turbulentScalarVelocityScale)
do iCell = 1, nCells

Expand Down Expand Up @@ -552,22 +555,29 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim
do i = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(i, iCell)

deltaVelocitySquared(1:minLevelCell(iCell)) = 0.0_RKIND
deltaVelocitySquared(1:minLevelEdgeBot(iEdge)) = 0.0_RKIND

normalVelocitySum(minLevelCell(iCell)) = normalVelocity(minLevelCell(iCell), iEdge)* &
layerThickEdgeMean(minLevelCell(iCell),iEdge)
layerThicknessEdgeSum(minLevelCell(iCell)) = layerThickEdgeMean(minLevelCell(iCell),iEdge)
normalVelocitySum(minLevelEdgeBot(iEdge)) = normalVelocity(minLevelEdgeBot(iEdge), iEdge)* &
layerThickEdgeMean(minLevelEdgeBot(iEdge),iEdge)
tangentialVelocitySum(minLevelEdgeBot(iEdge)) = tangentialVelocity(minLevelEdgeBot(iEdge),iEdge)* &
layerThickEdgeMean(minLevelEdgeBot(iEdge),iEdge)
layerThicknessEdgeSum(minLevelEdgeBot(iEdge)) = layerThickEdgeMean(minLevelEdgeBot(iEdge),iEdge)

do kIndexOBL = minLevelCell(iCell)+1, maxLevelCell(iCell)
do kIndexOBL = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge)
normalVelocitySum(kIndexOBL) = normalVelocitySum(kIndexOBL-1) + &
layerThickEdgeMean(kIndexOBL, iEdge)*normalVelocity(kIndexOBL, iEdge)
tangentialVelocitySum(kIndexOBL) = tangentialVelocitySum(kIndexOBL-1) + &
layerThickEdgeMean(kIndexOBL, iEdge)*tangentialVelocity(kIndexOBL, iEdge)
layerThicknessEdgeSum(kIndexOBL) = layerThicknessEdgeSum(kIndexOBL-1) + layerThickEdgeMean(kIndexOBL, iEdge)
end do

do kIndexOBL = minLevelCell(iCell), maxLevelCell(iCell)
do kIndexOBL = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
normalVelocityAv = normalVelocitySum(surfaceAverageIndex(kIndexOBL )) / &
layerThicknessEdgeSum(surfaceAverageIndex(kIndexOBL))
delU2 = ( normalVelocityAv - normalVelocity(kIndexOBL, iEdge) )**2
tangentialVelocityAv = tangentialVelocitySum(surfaceAverageIndex(kIndexOBL)) / &
layerThicknessEdgeSum(surfaceAverageIndex(kIndexOBL))
delU2 = (normalVelocityAv - normalVelocity(kIndexOBL, iEdge))**2 + &
(tangentialVelocityAv - tangentialVelocity(kIndexOBL, iEdge))**2
deltaVelocitySquared(kIndexOBL) = deltaVelocitySquared(kIndexOBL) + edgeAreaFractionOfCell(i,iCell) * delU2
end do
end do
Expand All @@ -584,7 +594,7 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim
do kIndexOBL = minLevelCell(iCell), maxLevelCell(iCell)
! !compute shear contribution assuming BLdepth is cell bottom
!Note that the factor of two is from averaging dot products to cell centers on a C-grid
bulkRichardsonNumberShear(kIndexOBL,iCell) = max(2.0_RKIND*deltaVelocitySquared(kIndexOBL), 1.0e-15_RKIND)
bulkRichardsonNumberShear(kIndexOBL,iCell) = max(deltaVelocitySquared(kIndexOBL), 1.0e-15_RKIND)

bulkRichardsonNumberBuoy(kIndexOBL,iCell) = gravity * (potentialDensity(kIndexOBL, iCell) &
- potentialDensitySum(surfaceAverageIndex(kIndexOBL)) &
Expand Down Expand Up @@ -832,6 +842,7 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim
deallocate(RiTemp)
deallocate(BVFSmoothed)
deallocate(normalVelocitySum)
deallocate(tangentialVelocitySum)
deallocate(potentialDensitySum)
deallocate(surfaceAverageIndex)
deallocate(deltaVelocitySquared)
Expand Down