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

Changes to richardson number calculation #66

Closed
wants to merge 6 commits into from
Closed

Changes to richardson number calculation #66

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
ed0432c
Changes to richardson number calculation
vanroekel Sep 6, 2023
75aa812
Fixes halo on tangentialVelocity
vanroekel Sep 15, 2023
4535ab7
Fixes to cvmix and domain decomp
vanroekel Sep 19, 2023
d9e3a22
Rolls back some extra halos
vanroekel Sep 19, 2023
486b05c
Adds halo exchange for RK4 and SI
vanroekel Sep 19, 2023
847acca
add domain to RK4 subroutine
mark-petersen Sep 19, 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
45 changes: 32 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 @@ -1595,29 +1597,46 @@ subroutine ocn_diagnostic_solve_richardson(displacedDensity, density, zMid, &
!$acc private(k, i, iEdge, delU2, shearSquared, shearMean)
#else
!$omp parallel
!$omp do schedule(runtime) private(k, shearSquared, i, iEdge, delU2, shearMean)
!$omp do schedule(runtime) private(k, shearSquared, i, iEdge, delU2, shearMean, &
!$omp jCell, shear, dz_edge, Ri_edge)
#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
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)
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 Expand Up @@ -2262,7 +2281,7 @@ subroutine ocn_diagnostic_solve_vortVel( &

! Compute tangential velocity

nEdges = nEdgesHalo(2)
vanroekel marked this conversation as resolved.
Show resolved Hide resolved
nEdges = nEdgesAll

#ifdef MPAS_OPENACC
!$acc parallel loop gang vector &
Expand Down
22 changes: 13 additions & 9 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 @@ -552,22 +552,26 @@ 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) &
+ 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)
layerThickEdgeMean(kIndexOBL, iEdge)*(normalVelocity(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
delU2 = ( normalVelocityAv - (tangentialVelocity(kIndexOBL, iEdge) + &
normalVelocity(kIndexOBL, iEdge)) )**2
deltaVelocitySquared(kIndexOBL) = deltaVelocitySquared(kIndexOBL) + edgeAreaFractionOfCell(i,iCell) * delU2
end do
end do
Expand All @@ -584,7 +588,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