Merge branch 'cbegeman/ocn/modify-zstar-for-inactive-top-cells' into …

…next (PR #5254)

Modify ocean z-star ALE coordinate for inactive top cells

This PR increases the generalizability of vertical coordinates by
relaxing the assumption that restingThickness spans ssh=0 to
bottomDepth. As such, it is tangentially related to inactive top cells.
The change affects cases where
  config_ALE_thickness_proportionality='restingThickness_times_weights'
which is the default.

[non-BFB]

components/mpas-ocean/src/shared/mpas_ocn_thick_ale.F

@@ -125,10 +125,12 @@ subroutine ocn_ALE_thickness(verticalMeshPool, SSH, ALE_thickness, &
          nCells        ! number of cells
 
       real (kind=RKIND) :: &
-         weightSum,    &! sum of weights in vertical
-         thicknessSum, &! total thickness
-         remainder,    &! track remainder in mix/max alteration
-         newThickness   ! temp used during min/max adjustment
+         weightSum,       &! sum of weights in vertical
+         restingThicknessSum,  &! total resting thickness
+         restingThicknessDiff, &! difference from total resting thickness
+         weightedThicknessSum, &! total resting thickness, weighted
+         remainder,       &! track remainder in mix/max alteration
+         newThickness      ! temp used during min/max adjustment
 
       real (kind=RKIND), dimension(:), allocatable :: &
          prelim_ALE_thickness,   & ! ALE thickness at new time
@@ -172,38 +174,44 @@ subroutine ocn_ALE_thickness(verticalMeshPool, SSH, ALE_thickness, &
          !xacc    present(ALE_thickness, SSH, restingThickness, &
          !xacc            minLevelCell, maxLevelCell, &
          !xacc            vertCoordMovementWeights) &
-         !xacc    private(k, kMin, kMax, thicknessSum)
+         !xacc    private(k, kMin, kMax, restingThicknessDiff, &
+         !xacc            restingThicknessSum, weightedThicknessSum)
 #else
          !$omp parallel
          !$omp do schedule(runtime) &
-         !$omp    private(k, kMin, kMax, thicknessSum)
+         !$omp    private(k, kMin, kMax, restingThicknessDiff, &
+         !$omp            restingThicknessSum, weightedThicknessSum)
 #endif
          do iCell = 1, nCells
             kMax = maxLevelCell(iCell)
             kMin = minLevelCell(iCell)
 
-            thicknessSum = 1e-14_RKIND
+            restingThicknessSum = 0.0_RKIND
+            weightedThicknessSum = 0.0_RKIND
             do k = kMin, kMax
-               thicknessSum = thicknessSum &
+               restingThicknessSum = restingThicknessSum &
+                            + restingThickness(k,iCell)
+               weightedThicknessSum = weightedThicknessSum &
                             + vertCoordMovementWeights(k) &
                             * restingThickness(k,iCell)
             end do
 
+            restingThicknessDiff = restingThicknessSum - bottomDepth(iCell)
             ! Note that restingThickness is nonzero, and remaining
             ! terms are perturbations about zero.
             ! This is equation 4 and 6 in Petersen et al 2015,
             ! but with eqn 6
             do k = kMin, kMax
-               ALE_thickness(k,iCell) = restingThickness(k,iCell) &
-                  + (SSH(iCell)*vertCoordMovementWeights(k)* &
-                     restingThickness(k,iCell) )/thicknessSum
+               ALE_thickness(k, iCell) = restingThickness(k, iCell) &
+                  + ( (SSH(iCell) + restingThicknessDiff) * vertCoordMovementWeights(k) * &
+                      restingThickness(k, iCell) ) / max(1e-14_RKIND, weightedThicknessSum)
             end do
          enddo
 #ifndef MPAS_OPENACC
          !$omp end do
          !$omp end parallel
 #endif
-   
+
       ! weights_only
       case (ALEthickProportionWgtsOnly)