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 417ed8f91014..8a3c633536c0 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 @@ -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 @@ -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 @@ -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) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F index 5e145dff3f32..743c250ed306 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F @@ -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) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index 0fd6f1102533..ef259b575eba 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -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) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index 5d3df864d25e..f9368e035719 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -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}) ! @@ -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 - ! 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!}}} !*********************************************************************** diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vmix_cvmix.F b/components/mpas-ocean/src/shared/mpas_ocn_vmix_cvmix.F index 564d8ded9311..fcede65e420c 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vmix_cvmix.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vmix_cvmix.F @@ -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 @@ -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)) @@ -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 @@ -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 @@ -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)) & @@ -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)