Skip to content

Commit

Permalink
Merge branch 'vanroekel/ocean/fix-richardson-number-calculation' (PR #…
Browse files Browse the repository at this point in the history 
…5946)

Fixes MPAS-Ocean Richardson number calculation

This changes the calculation of the gradient and bulk richardson numbers
to utilize the normal and tangential velocities to avoid the factor of
two issue and also fixes a number of places where the loop was over
maxLevelCell and should have been to maxLevelEdgeTop.

[NML]
[CC]
  • Loading branch information
@jonbob
jonbob committed Oct 26, 2023
2 parents 4d48d07 + cfb1208 commit 510fca6
Show file tree
Hide file tree
Showing 5 changed files with 63 additions and 28 deletions.
7 changes: 5 additions & 2 deletions components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F
View file
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
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
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
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
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

0 comments on commit 510fca6

Please sign in to comment.