Fix bounds on DivHU computation to avoid inactive layers

Author: xylar

components/omega/src/ocn/VertAdv.cpp

@@ -379,6 +379,7 @@ void VertAdv::computeVerticalPseudoVelocity(
    OMEGA_SCOPE(LocEOnC, Mesh->EdgesOnCell);
    OMEGA_SCOPE(LocDvE, Mesh->DvEdge);
    OMEGA_SCOPE(LocESOnC, Mesh->EdgeSignOnCell);
+   OMEGA_SCOPE(LocMaxLayerEdgeTop, VCoord->MaxLayerEdgeTop);
 
    // Loop over all cells owned by the task
    parallelForOuter(
@@ -402,12 +403,16 @@ void VertAdv::computeVerticalPseudoVelocity(
                  const I4 KLen            = chunkLength(KChunk, KStart, KMax);
 
                  for (int J = 0; J < LocNEOnC(ICell); ++J) {
-                    const I4 JEdge = LocEOnC(ICell, J);
+                    const I4 JEdge    = LocEOnC(ICell, J);
+                    const I4 MaxKEdge = LocMaxLayerEdgeTop(JEdge);
                     for (int KVec = 0; KVec < KLen; ++KVec) {
                        const I4 K = KStart + KVec;
-                       DivHUTmp[KVec] -= LocDvE(JEdge) * LocESOnC(ICell, J) *
-                                         FluxPseudoThickEdge(JEdge, K) *
-                                         NormalVelocity(JEdge, K) * InvAreaCell;
+                       if (K <= MaxKEdge) {
+                          DivHUTmp[KVec] -= LocDvE(JEdge) * LocESOnC(ICell, J) *
+                                            FluxPseudoThickEdge(JEdge, K) *
+                                            NormalVelocity(JEdge, K) *
+                                            InvAreaCell;
+                       }
                     }
                  }
                  for (int KVec = 0; KVec < KLen; ++KVec) {