From f9d25732a10ec8ce86f922fb850e5f7b401606be Mon Sep 17 00:00:00 2001
From: Carolyn Begeman <cbegeman@lanl.gov>
Date: Wed, 29 Mar 2023 09:28:51 -0500
Subject: [PATCH] Use alternate weighted z-star formulation

---
 .../src/shared/mpas_ocn_thick_ale.F           | 39 ++++++++++++-------
 1 file changed, 25 insertions(+), 14 deletions(-)

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_thick_ale.F b/components/mpas-ocean/src/shared/mpas_ocn_thick_ale.F
index a4e861f6ed67..a99737b4b755 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_thick_ale.F
+++ b/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 over 1 to maxLevelCell
+         activeWeightSum,    &! sum of weights in vertical over minLevelCell to maxLevelCell
+         restingThicknessSum,&! total resting thickness
+         newThicknessSum,    &! total thickness at current time
+         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
@@ -169,34 +171,43 @@ subroutine ocn_ALE_thickness(verticalMeshPool, SSH, ALE_thickness, &
 
 #ifdef MPAS_OPENACC
          !xacc parallel loop &
-         !xacc    present(ALE_thickness, SSH, restingThickness, &
+         !xacc    present(ALE_thickness, SSH, restingThicknessSum, &
+         !xacc            newThicknessSum, weightSum, activeWeightSum, &
          !xacc            minLevelCell, maxLevelCell, &
          !xacc            vertCoordMovementWeights) &
-         !xacc    private(k, kMin, kMax, thicknessSum)
+         !xacc    private(k, kMin, kMax, restingThicknessSum, &
+         !xacc            newThicknessSum, weightSum, activeWeightSum)
 #else
          !$omp parallel
          !$omp do schedule(runtime) &
-         !$omp    private(k, kMin, kMax, thicknessSum)
+         !$omp    private(k, kMin, kMax, restingThicknessSum, &
+         !$omp            newThicknessSum, weightSum, activeWeightSum)
 #endif
          do iCell = 1, nCells
             kMax = maxLevelCell(iCell)
             kMin = minLevelCell(iCell)
    
-            thicknessSum = 1e-14_RKIND
+            restingThicknessSum = 1e-14_RKIND
+            weightSum = 0.0_RKIND
+            activeWeightSum = 1e-14_RKIND
             do k = kMin, kMax
-               thicknessSum = thicknessSum &
-                            + vertCoordMovementWeights(k) &
-                            * restingThickness(k,iCell)
+               restingThicknessSum = restingThicknessSum &
+                            + restingThickness(k,iCell)
+               activeWeightSum = activeWeightSum + vertCoordMovementWeights(k)
             end do
+            do k = 1, kMax
+               weightSum = weightSum + vertCoordMovementWeights(k)
+            end do
+            newThicknessSum = bottomDepth(iCell) + ssh(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) = (newThicknessSum / restingThicknessSum ) &
+                      * (weightSum / activeWeightSum) &
+                      * vertCoordMovementWeights(k) * restingThickness(k, iCell)
             end do
          enddo
 #ifndef MPAS_OPENACC