Better variables for updated masks in fill_miss_2d

  Within fill_miss_2d, there are 2 sets of masked values used for remapping -
masked values in neighboring points in arbitrary units (with names like "south")
and the updated valid mask with values of 0 or 1 (with names like "gs").  At one
point in the code, the former were being used where the latter were more
appropriate.  This commit changes this to help with the understandability of the
code.

  However, when the order of arithmetic is not specified with parentheses, the
Nvidia compiler is changing answers with this change of variables.  To avoid
this change in answers, this commit is only changing which variables are used
when the answer_date is larger than 20190101 and the order of the sums are
specified with parentheses.  This commit also now moves the logical test for
ans_2018 outside of the i- and j- loops where non-zero values of a_chg are set.

  All answers and output are bitwise identical.

src/framework/MOM_horizontal_regridding.F90

@@ -125,7 +125,7 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, acrit, num_pass, relc, debug,
   real, dimension(SZI_(G),SZJ_(G)) :: good_new ! The values of good_ to use for the next iteration [nondim]
 
   real    :: east, west, north, south ! Valid neighboring values or 0 for invalid values [arbitrary]
-  real    :: ge, gw, gn, gs  ! Flags indicating which neighbors have valid values [nondim]
+  real    :: ge, gw, gn, gs  ! Flags set to 0 or 1 indicating which neighbors have valid values [nondim]
   real    :: ngood     ! The number of valid values in neighboring points [nondim]
   real    :: nfill     ! The remaining number of points to fill [nondim]
   real    :: nfill_prev ! The previous value of nfill [nondim]
@@ -227,23 +227,30 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, acrit, num_pass, relc, debug,
   ! Do Laplacian smoothing for the points that have been filled in.
   do k=1,npass
     call pass_var(aout,G%Domain)
-    do j=js,je ; do i=is,ie
-      if (fill(i,j) == 1) then
-        east = max(good(i+1,j),fill(i+1,j)) ; west = max(good(i-1,j),fill(i-1,j))
-        north = max(good(i,j+1),fill(i,j+1)) ; south = max(good(i,j-1),fill(i,j-1))
-        if (ans_2018) then
+
+    a_chg(:,:) = 0.0
+    if (ans_2018) then
+      do j=js,je ; do i=is,ie
+        if (fill(i,j) == 1) then
+          east = max(good(i+1,j),fill(i+1,j)) ; west = max(good(i-1,j),fill(i-1,j))
+          north = max(good(i,j+1),fill(i,j+1)) ; south = max(good(i,j-1),fill(i,j-1))
           a_chg(i,j) = relax_coeff*(south*aout(i,j-1)+north*aout(i,j+1) + &
                                     west*aout(i-1,j)+east*aout(i+1,j) - &
                                    (south+north+west+east)*aout(i,j))
-        else
-          a_chg(i,j) = relax_coeff*( ((south*aout(i,j-1) + north*aout(i,j+1)) + &
-                                  (west*aout(i-1,j)+east*aout(i+1,j))) - &
-                                 ((south+north)+(west+east))*aout(i,j) )
         endif
-      else
-        a_chg(i,j) = 0.
-      endif
-    enddo ; enddo
+      enddo ; enddo
+    else
+      do j=js,je ; do i=is,ie
+        if (fill(i,j) == 1) then
+          ge = max(good(i+1,j),fill(i+1,j)) ; gw = max(good(i-1,j),fill(i-1,j))
+          gn = max(good(i,j+1),fill(i,j+1)) ; gs = max(good(i,j-1),fill(i,j-1))
+          a_chg(i,j) = relax_coeff*( ((gs*aout(i,j-1) + gn*aout(i,j+1)) + &
+                                      (gw*aout(i-1,j) + ge*aout(i+1,j))) - &
+                                     ((gs + gn) + (gw + ge))*aout(i,j) )
+        endif
+      enddo ; enddo
+    endif
+
     ares = 0.0
     do j=js,je ; do i=is,ie
       aout(i,j) = a_chg(i,j) + aout(i,j)