Merge branch 'tcclevenger/eamxx/fpe_shoc_fix' (PR E3SM-Project#7368)

Fixes FPE introduced in E3SM-Project#7188 from temporary values being computed
which attempt to take sqrt(some negative value). These values weren't
being used, so not an actual bug.

[BFB]

Author: tcclevenger

components/eamxx/src/physics/shoc/impl/shoc_compute_shoc_mix_shoc_length_impl.hpp

@@ -25,38 +25,39 @@ ::compute_shoc_mix_shoc_length(
   const Int nlev_pack = ekat::npack<Spack>(nlev);
   const auto maxlen = scream::shoc::Constants<Scalar>::maxlen;
   const auto vk = C::Karman;
-  
+
   // Eddy turnover timescale
   const Scalar tscale = 400;
 
   Kokkos::parallel_for(Kokkos::TeamVectorRange(team, nlev_pack), [&] (const Int& k) {
     const Spack tkes = ekat::sqrt(tke(k));
     const Spack brunt2 = ekat::max(0, brunt(k));
 
-   if (shoc_1p5tke){
-
-        // If 1.5 TKE closure then set length scale to vertical grid spacing for
-	//   cells with unstable brunt vaisalla frequency.  Otherwise, overwrite the length
-	//   scale in stable cells with the new definition.
+    if (shoc_1p5tke){
+      // If 1.5 TKE closure then set length scale to vertical grid spacing for
+	    //   cells with unstable brunt vaisalla frequency.  Otherwise, overwrite the length
+	    //   scale in stable cells with the new definition.
 
-	// Search for stable cells
-        const auto stable_mask = brunt(k) > 0;
+	    // Search for stable cells
+      const auto stable_mask = brunt(k) > 0;
 
-        // Define length scale for stable cells
-        const auto length_tmp = ekat::sqrt(0.76*tk(k)/0.1/ekat::sqrt(brunt(k) + 1.e-10));
-        // Limit the stability corrected length scale between 0.1*dz and dz
-	const auto limited_len = ekat::min(dz_zt(k),ekat::max(0.1*dz_zt(k),length_tmp));
+      // To avoid FPE when calculating sqrt(brunt), set brunt_tmp=0 in the case brunt<1.
+      Spack brunt_tmp(stable_mask, brunt(k));
 
-        // Set length scale to vertical grid if unstable, otherwise the stability adjusted value.
-        shoc_mix(k).set(stable_mask, limited_len, dz_zt(k));
+      // Define length scale for stable cells
+      const auto length_tmp = ekat::sqrt(0.76*tk(k)/0.1/ekat::sqrt(brunt_tmp + 1.e-10));
+      // Limit the stability corrected length scale between 0.1*dz and dz
+	    const auto limited_len = ekat::min(dz_zt(k),ekat::max(0.1*dz_zt(k),length_tmp));
 
-     }else{
-        shoc_mix(k) = ekat::min(maxlen,
-                            sp(2.8284)*(ekat::sqrt(1/((1/(tscale*tkes*vk*zt_grid(k)))
-                            + (1/(tscale*tkes*l_inf))
-                            + sp(0.01)*(brunt2/tke(k)))))/length_fac);
+      // Set length scale to vertical grid if unstable, otherwise the stability adjusted value.
+      shoc_mix(k).set(stable_mask, limited_len, dz_zt(k));
+    } else{
+      shoc_mix(k) = ekat::min(maxlen,
+                              sp(2.8284)*(ekat::sqrt(1/((1/(tscale*tkes*vk*zt_grid(k)))
+                              + (1/(tscale*tkes*l_inf))
+                              + sp(0.01)*(brunt2/tke(k)))))/length_fac);
 
-     }
+    }
   });
 }