Merge branch 'vanroekel/ocean/fix-richardson-number-calculation' into…

… next (PR #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. [CC]
E3SM-Project · Oct 24, 2023 · 444b05f · 444b05f
2 parents 5472a90 + 5c6f47c
commit 444b05f
Show file tree

Hide file tree

Showing 5 changed files with 61 additions and 27 deletions.
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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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)