Implement changes suggested by @Hallberg-NOAA

On PR#670, @Hallberg-NOAA pointed out two aspects of the code that could be
changed. As detailed in a reply to those comments, one of these was a
non-answer changing bit of logic that could safely be removed and the other
involved changing array syntax arithmetic to explicit loops.

src/tracer/MOM_neutral_diffusion.F90

@@ -389,6 +389,9 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS)
   integer :: iMethod
   real, dimension(SZI_(G)) :: ref_pres ! Reference pressure used to calculate alpha/beta
   real :: h_neglect, h_neglect_edge
+  real :: pa_to_H
+
+  pa_to_H = 1. / GV%H_to_pa
 
   !### Try replacing both of these with GV%H_subroundoff
   if (GV%Boussinesq) then
@@ -444,15 +447,13 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS)
     else ! Discontinuous reconstruction
       do k = 1, G%ke
         if (CS%ref_pres<0) ref_pres(:) = CS%Pint(:,j,k)
-        if (CS%stable_cell(i,j,k)) &
-          call calculate_density_derivs(CS%T_i(:,j,k,1), CS%S_i(:,j,k,1), ref_pres, &
-                                        CS%dRdT_i(:,j,k,1), CS%dRdS_i(:,j,k,1), G%isc-1, G%iec-G%isc+3, CS%EOS)
+        call calculate_density_derivs(CS%T_i(:,j,k,1), CS%S_i(:,j,k,1), ref_pres, &
+                                      CS%dRdT_i(:,j,k,1), CS%dRdS_i(:,j,k,1), G%isc-1, G%iec-G%isc+3, CS%EOS)
         if (CS%ref_pres<0) then
           ref_pres(:) = CS%Pint(:,j,k+1)
         endif
-        if (CS%stable_cell(i,j,k)) &
-          call calculate_density_derivs(CS%T_i(:,j,k,2), CS%S_i(:,j,k,2), ref_pres, &
-                                        CS%dRdT_i(:,j,k,2), CS%dRdS_i(:,j,k,2), G%isc-1, G%iec-G%isc+3, CS%EOS)
+        call calculate_density_derivs(CS%T_i(:,j,k,2), CS%S_i(:,j,k,2), ref_pres, &
+                                      CS%dRdT_i(:,j,k,2), CS%dRdS_i(:,j,k,2), G%isc-1, G%iec-G%isc+3, CS%EOS)
       enddo
     endif
   enddo
@@ -518,9 +519,16 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS)
     endif
   enddo ; enddo
 
+  ! Continuous reconstructions calculate hEff as the difference between the pressures of the neutral surfaces which
+  ! need to be reconverted to thickness units. The discontinuous version calculates hEff from the fraction of the
+  ! nondimensional fraction of the layer occupied by the
   if (CS%continuous_reconstruction) then
-    CS%uhEff(:,:,:) = CS%uhEff(:,:,:) / GV%H_to_pa
-    CS%vhEff(:,:,:) = CS%vhEff(:,:,:) / GV%H_to_pa
+    do k = 1, CS%nsurf-1 ; do j = G%jsc, G%jec ; do I = G%isc-1, G%iec
+      if (G%mask2dCu(I,j) > 0.) CS%uhEff(I,j,k) = CS%uhEff(I,j,k) * pa_to_H
+    enddo ; enddo ; enddo
+    do k = 1, CS%nsurf-1 ; do J = G%jsc-1, G%jec ; do i = G%isc, G%iec
+      if (G%mask2dCv(i,J) > 0.) CS%vhEff(i,J,k) = CS%vhEff(i,J,k) * pa_to_H
+    enddo ; enddo ; enddo
   endif
 
   if (CS%id_uhEff_2d>0) then

src/tracer/MOM_neutral_diffusion_aux.F90

@@ -161,7 +161,7 @@ real function calc_drho(T1, S1, dRdT1, dRdS1, T2, S2, dRdT2, dRdS2)
 end function calc_drho
 
 !> Calculate the difference in neutral density between a reference T, S, alpha, and beta
-!! and a poiet on the polynomial reconstructions of T, S
+!! at a point on the polynomial reconstructions of T, S
 subroutine drho_at_pos(CS, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, x0, &
                        delta_rho, P_out, T_out, S_out, alpha_avg_out, beta_avg_out, delta_T_out, delta_S_out)
   type(ndiff_aux_CS_type), intent(in) :: CS           !< Control structure with parameters for this module