From 3f041d93fbb784b420456ca2e6b60df647133426 Mon Sep 17 00:00:00 2001
From: Marshall Ward <marshall.ward@noaa.gov>
Date: Thu, 3 Oct 2019 18:24:14 -0400
Subject: [PATCH] MEKE diagnostic array fixes

This patch fixes the following MEKE diagnostics:

- MEKE_Ue, MEKE_Ub, MEKE_Ut

The diagnostics were computed as inline operations inside post_data,
e.g.:

    post_data(..., sqrt(0, max(0., MEKE*bottomFac2)))

rather than computing the fields explicitly inside of array loops.

This case causing floating point exceptions in Intel compilers, possibly
likely due to evaluations inside of halos.

We resolve these diagnostics by computing the values into a scratch
array which is then passed to post_data.
---
 src/parameterizations/lateral/MOM_MEKE.F90 | 28 ++++++++++++++++++----
 1 file changed, 24 insertions(+), 4 deletions(-)

diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90
index 853f3a8613..8d3e78262c 100644
--- a/src/parameterizations/lateral/MOM_MEKE.F90
+++ b/src/parameterizations/lateral/MOM_MEKE.F90
@@ -139,7 +139,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
     del4MEKE, &     ! Time-integrated MEKE tendency arising from the biharmonic of MEKE [L2 T-2 ~> m2 s-2].
     LmixScale, &    ! Eddy mixing length [L ~> m].
     barotrFac2, &   ! Ratio of EKE_barotropic / EKE [nondim]
-    bottomFac2      ! Ratio of EKE_bottom / EKE [nondim]
+    bottomFac2, &   ! Ratio of EKE_bottom / EKE [nondim]
+    tmp             ! Temporary variable for diagnostic computation
 
   real, dimension(SZIB_(G),SZJ_(G)) :: &
     MEKE_uflux, &   ! The zonal advective and diffusive flux of MEKE with different units in different
@@ -593,10 +594,29 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
     endif
 
     ! Offer fields for averaging.
+
+    if (any([CS%id_Ue, CS%id_Ub, CS%id_Ut] > 0)) &
+      tmp(:,:) = 0.
+
     if (CS%id_MEKE>0) call post_data(CS%id_MEKE, MEKE%MEKE, CS%diag)
-    if (CS%id_Ue>0) call post_data(CS%id_Ue, sqrt(max(0.,2.0*MEKE%MEKE)), CS%diag)
-    if (CS%id_Ub>0) call post_data(CS%id_Ub, sqrt(max(0.,2.0*MEKE%MEKE*bottomFac2)), CS%diag)
-    if (CS%id_Ut>0) call post_data(CS%id_Ut, sqrt(max(0.,2.0*MEKE%MEKE*barotrFac2)), CS%diag)
+    if (CS%id_Ue>0) then
+      do j=js,je ; do i=is,ie
+        tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j)))
+      enddo ; enddo
+      call post_data(CS%id_Ue, tmp, CS%diag)
+    endif
+    if (CS%id_Ub>0) then
+      do j=js,je ; do i=is,ie
+        tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j) * bottomFac2(i,j)))
+      enddo ; enddo
+      call post_data(CS%id_Ub, tmp, CS%diag)
+    endif
+    if (CS%id_Ut>0) then
+      do j=js,je ; do i=is,ie
+        tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j) * barotrFac2(i,j)))
+      enddo ; enddo
+      call post_data(CS%id_Ut, tmp, CS%diag)
+    endif
     if (CS%id_Kh>0) call post_data(CS%id_Kh, MEKE%Kh, CS%diag)
     if (CS%id_Ku>0) call post_data(CS%id_Ku, MEKE%Ku, CS%diag)
     if (CS%id_Au>0) call post_data(CS%id_Au, MEKE%Au, CS%diag)