diff --git a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
index acbbc292de..faa74a7fe0 100644
--- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
+++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90
@@ -374,7 +374,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
       fluxes%salt_flux_added(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) ! Diagnostic
     else
       do j=js,je ; do i=is,ie
-        if (G%mask2dT(i,j) > 0.5) then
+        if (G%mask2dT(i,j) > 0.0) then
           delta_sss = sfc_state%SSS(i,j) - data_restore(i,j)
           delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),CS%max_delta_srestore)
           fluxes%vprec(i,j) = (CS%basin_mask(i,j)*open_ocn_mask(i,j)*CS%srestore_mask(i,j))* &
@@ -956,13 +956,13 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
       if (present(taux).and.present(tauy)) then
         do j=jsh,jeh ; do I=Isqh,Ieqh
           taux(I,j) = 0.0
-          if ((G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) > 0) &
+          if ((G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) > 0.0) &
             taux(I,j) = (G%mask2dBu(I,J)*taux_in_B(I,J) + G%mask2dBu(I,J-1)*taux_in_B(I,J-1)) / &
                         (G%mask2dBu(I,J) + G%mask2dBu(I,J-1))
         enddo ; enddo
         do J=Jsqh,Jeqh ; do i=ish,ieh
           tauy(i,J) = 0.0
-          if ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) > 0) &
+          if ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) > 0.0) &
             tauy(i,J) = (G%mask2dBu(I,J)*tauy_in_B(I,J) + G%mask2dBu(I-1,J)*tauy_in_B(I-1,J)) / &
                         (G%mask2dBu(I,J) + G%mask2dBu(I-1,J))
         enddo ; enddo
@@ -984,14 +984,14 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
 
       if (present(taux)) then ; do j=jsh,jeh ; do I=Isqh,Ieqh
         taux(I,j) = 0.0
-        if ((G%mask2dT(i,j) + G%mask2dT(i+1,j)) > 0) &
+        if ((G%mask2dT(i,j) + G%mask2dT(i+1,j)) > 0.0) &
           taux(I,j) = (G%mask2dT(i,j)*taux_in_A(i,j) + G%mask2dT(i+1,j)*taux_in_A(i+1,j)) / &
                       (G%mask2dT(i,j) + G%mask2dT(i+1,j))
       enddo ; enddo ; endif
 
       if (present(tauy)) then ; do J=Jsqh,Jeqh ; do i=ish,ieh
         tauy(i,J) = 0.0
-        if ((G%mask2dT(i,j) + G%mask2dT(i,j+1)) > 0) &
+        if ((G%mask2dT(i,j) + G%mask2dT(i,j+1)) > 0.0) &
           tauy(i,J) = (G%mask2dT(i,j)*tauy_in_A(i,j) + G%mask2dT(i,J+1)*tauy_in_A(i,j+1)) / &
                       (G%mask2dT(i,j) + G%mask2dT(i,j+1))
       enddo ; enddo ; endif
@@ -1029,10 +1029,10 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
         gustiness = CS%gust_const
         if (CS%read_gust_2d) then
           if ((wind_stagger == CGRID_NE) .or. &
-              ((wind_stagger == AGRID) .and. (G%mask2dT(i,j) > 0)) .or. &
+              ((wind_stagger == AGRID) .and. (G%mask2dT(i,j) > 0.0)) .or. &
               ((wind_stagger == BGRID_NE) .and. &
                (((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + &
-                (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0)) ) &
+                (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0.0)) ) &
             gustiness = CS%gust(i,j)
         endif
         ustar(i,j) = sqrt(gustiness*IRho0 + IRho0*Pa_conversion*IOB%stress_mag(i-i0,j-j0))
@@ -1050,7 +1050,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
       do j=js,je ; do i=is,ie
         tau_mag = 0.0 ; gustiness = CS%gust_const
         if (((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + &
-             (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0) then
+             (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0.0) then
           tau_mag = sqrt(((G%mask2dBu(I,J)*(taux_in_B(I,J)**2 + tauy_in_B(I,J)**2) + &
               G%mask2dBu(I-1,J-1)*(taux_in_B(I-1,J-1)**2 + tauy_in_B(I-1,J-1)**2)) + &
              (G%mask2dBu(I,J-1)*(taux_in_B(I,J-1)**2 + tauy_in_B(I,J-1)**2) + &
@@ -1069,7 +1069,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
       do j=js,je ; do i=is,ie
         tau_mag = G%mask2dT(i,j) * sqrt(taux_in_A(i,j)**2 + tauy_in_A(i,j)**2)
         gustiness = CS%gust_const
-        if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0)) gustiness = CS%gust(i,j)
+        if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0.0)) gustiness = CS%gust(i,j)
         if (do_ustar) ustar(i,j) = sqrt(gustiness*IRho0 + IRho0 * tau_mag)
         if (CS%answers_2018) then
           if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0)
@@ -1080,10 +1080,10 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
     else  ! C-grid wind stresses.
       do j=js,je ; do i=is,ie
         taux2 = 0.0 ; tauy2 = 0.0
-        if ((G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) > 0) &
+        if ((G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) > 0.0) &
           taux2 = (G%mask2dCu(I-1,j)*taux_in_C(I-1,j)**2 + G%mask2dCu(I,j)*taux_in_C(I,j)**2) / &
                   (G%mask2dCu(I-1,j) + G%mask2dCu(I,j))
-        if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0) &
+        if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0.0) &
           tauy2 = (G%mask2dCv(i,J-1)*tauy_in_C(i,J-1)**2 + G%mask2dCv(i,J)*tauy_in_C(i,J)**2) / &
                   (G%mask2dCv(i,J-1) + G%mask2dCv(i,J))
         tau_mag = sqrt(taux2 + tauy2)
diff --git a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90
index 7b32270b4c..a8a3ea0c69 100644
--- a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90
+++ b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90
@@ -385,7 +385,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
       fluxes%salt_flux_added(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) ! Diagnostic
     else
       do j=js,je ; do i=is,ie
-        if (G%mask2dT(i,j) > 0.5) then
+        if (G%mask2dT(i,j) > 0.0) then
           delta_sss = sfc_state%SSS(i,j) - data_restore(i,j)
           delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),CS%max_delta_srestore)
           fluxes%vprec(i,j) = (CS%basin_mask(i,j)*open_ocn_mask(i,j)*CS%srestore_mask(i,j))* &
@@ -749,7 +749,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     do j=js,je ; do I=Isq,Ieq
       forces%taux(I,j) = 0.0
-      if ((G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) > 0) &
+      if ((G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) > 0.0) &
         forces%taux(I,j) = (G%mask2dBu(I,J)*taux_at_q(I,J) + &
                             G%mask2dBu(I,J-1)*taux_at_q(I,J-1)) / &
                            (G%mask2dBu(I,J) + G%mask2dBu(I,J-1))
@@ -757,7 +757,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     do J=Jsq,Jeq ; do i=is,ie
       forces%tauy(i,J) = 0.0
-      if ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) > 0) &
+      if ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) > 0.0) &
         forces%tauy(i,J) = (G%mask2dBu(I,J)*tauy_at_q(I,J) + &
                             G%mask2dBu(I-1,J)*tauy_at_q(I-1,J)) / &
                            (G%mask2dBu(I,J) + G%mask2dBu(I-1,J))
@@ -770,7 +770,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
     do j=js,je ; do i=is,ie
       tau_mag = 0.0 ; gustiness = CS%gust_const
       if (((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + &
-           (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0) then
+           (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0.0) then
         tau_mag = sqrt(((G%mask2dBu(I,J)*(taux_at_q(I,J)**2 + tauy_at_q(I,J)**2) + &
              G%mask2dBu(I-1,J-1)*(taux_at_q(I-1,J-1)**2 + tauy_at_q(I-1,J-1)**2)) + &
              (G%mask2dBu(I,J-1)*(taux_at_q(I,J-1)**2 + tauy_at_q(I,J-1)**2) + &
@@ -786,7 +786,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     do j=js,je ; do I=Isq,Ieq
       forces%taux(I,j) = 0.0
-      if ((G%mask2dT(i,j) + G%mask2dT(i+1,j)) > 0) &
+      if ((G%mask2dT(i,j) + G%mask2dT(i+1,j)) > 0.0) &
         forces%taux(I,j) = (G%mask2dT(i,j)*taux_at_h(i,j) + &
                             G%mask2dT(i+1,j)*taux_at_h(i+1,j)) / &
                            (G%mask2dT(i,j) + G%mask2dT(i+1,j))
@@ -794,7 +794,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     do J=Jsq,Jeq ; do i=is,ie
       forces%tauy(i,J) = 0.0
-      if ((G%mask2dT(i,j) + G%mask2dT(i,j+1)) > 0) &
+      if ((G%mask2dT(i,j) + G%mask2dT(i,j+1)) > 0.0) &
         forces%tauy(i,J) = (G%mask2dT(i,j)*tauy_at_h(i,j) + &
                             G%mask2dT(i,J+1)*tauy_at_h(i,j+1)) / &
                            (G%mask2dT(i,j) + G%mask2dT(i,j+1))
@@ -802,7 +802,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     do j=js,je ; do i=is,ie
       gustiness = CS%gust_const
-      if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0)) gustiness = CS%gust(i,j)
+      if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0.0)) gustiness = CS%gust(i,j)
       forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * &
                                sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2))
     enddo; enddo
@@ -814,12 +814,12 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     do j=js,je ; do i=is,ie
       taux2 = 0.0
-      if ((G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) > 0) &
+      if ((G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) > 0.0) &
         taux2 = (G%mask2dCu(I-1,j)*forces%taux(I-1,j)**2 + &
                  G%mask2dCu(I,j)*forces%taux(I,j)**2) / (G%mask2dCu(I-1,j) + G%mask2dCu(I,j))
 
       tauy2 = 0.0
-      if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0) &
+      if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0.0) &
         tauy2 = (G%mask2dCv(i,J-1)*forces%tauy(i,J-1)**2 + &
                  G%mask2dCv(i,J)*forces%tauy(i,J)**2) / (G%mask2dCv(i,J-1) + G%mask2dCv(i,J))
 
diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
index c704214930..d04ecf10e4 100644
--- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
+++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
@@ -410,7 +410,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
       fluxes%salt_flux_added(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) ! Diagnostic
     else
       do j=js,je ; do i=is,ie
-        if (G%mask2dT(i,j) > 0.5) then
+        if (G%mask2dT(i,j) > 0.0) then
           delta_sss = sfc_state%SSS(i,j) - data_restore(i,j)
           delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),CS%max_delta_srestore)
           fluxes%vprec(i,j) = (CS%basin_mask(i,j)*open_ocn_mask(i,j)*CS%srestore_mask(i,j))* &
@@ -806,7 +806,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     do j=js,je ; do I=Isq,Ieq
       forces%taux(I,j) = 0.0
-      if ((G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) > 0) &
+      if ((G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) > 0.0) &
         forces%taux(I,j) = (G%mask2dBu(I,J)*taux_at_q(I,J) + &
                             G%mask2dBu(I,J-1)*taux_at_q(I,J-1)) / &
                            (G%mask2dBu(I,J) + G%mask2dBu(I,J-1))
@@ -814,7 +814,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     do J=Jsq,Jeq ; do i=is,ie
       forces%tauy(i,J) = 0.0
-      if ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) > 0) &
+      if ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) > 0.0) &
         forces%tauy(i,J) = (G%mask2dBu(I,J)*tauy_at_q(I,J) + &
                             G%mask2dBu(I-1,J)*tauy_at_q(I-1,J)) / &
                            (G%mask2dBu(I,J) + G%mask2dBu(I-1,J))
@@ -827,7 +827,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
     do j=js,je ; do i=is,ie
       tau_mag = 0.0 ; gustiness = CS%gust_const
       if (((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + &
-           (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0) then
+           (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0.0) then
         tau_mag = sqrt(((G%mask2dBu(I,J)*(taux_at_q(I,J)**2 + tauy_at_q(I,J)**2) + &
             G%mask2dBu(I-1,J-1)*(taux_at_q(I-1,J-1)**2 + tauy_at_q(I-1,J-1)**2)) + &
            (G%mask2dBu(I,J-1)*(taux_at_q(I,J-1)**2 + tauy_at_q(I,J-1)**2) + &
@@ -843,7 +843,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     do j=js,je ; do I=Isq,Ieq
       forces%taux(I,j) = 0.0
-      if ((G%mask2dT(i,j) + G%mask2dT(i+1,j)) > 0) &
+      if ((G%mask2dT(i,j) + G%mask2dT(i+1,j)) > 0.0) &
         forces%taux(I,j) = (G%mask2dT(i,j)*taux_at_h(i,j) + &
                             G%mask2dT(i+1,j)*taux_at_h(i+1,j)) / &
                            (G%mask2dT(i,j) + G%mask2dT(i+1,j))
@@ -851,7 +851,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     do J=Jsq,Jeq ; do i=is,ie
       forces%tauy(i,J) = 0.0
-      if ((G%mask2dT(i,j) + G%mask2dT(i,j+1)) > 0) &
+      if ((G%mask2dT(i,j) + G%mask2dT(i,j+1)) > 0.0) &
         forces%tauy(i,J) = (G%mask2dT(i,j)*tauy_at_h(i,j) + &
                             G%mask2dT(i,J+1)*tauy_at_h(i,j+1)) / &
                            (G%mask2dT(i,j) + G%mask2dT(i,j+1))
@@ -859,7 +859,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     do j=js,je ; do i=is,ie
       gustiness = CS%gust_const
-      if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0)) gustiness = CS%gust(i,j)
+      if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0.0)) gustiness = CS%gust(i,j)
       forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * &
                                sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2))
     enddo ; enddo
@@ -871,12 +871,12 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     do j=js,je ; do i=is,ie
       taux2 = 0.0
-      if ((G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) > 0) &
+      if ((G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) > 0.0) &
         taux2 = (G%mask2dCu(I-1,j)*forces%taux(I-1,j)**2 + &
                  G%mask2dCu(I,j)*forces%taux(I,j)**2) / (G%mask2dCu(I-1,j) + G%mask2dCu(I,j))
 
       tauy2 = 0.0
-      if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0) &
+      if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0.0) &
         tauy2 = (G%mask2dCv(i,J-1)*forces%tauy(i,J-1)**2 + &
                  G%mask2dCv(i,J)*forces%tauy(i,J)**2) / (G%mask2dCv(i,J-1) + G%mask2dCv(i,J))
 
diff --git a/config_src/drivers/solo_driver/MESO_surface_forcing.F90 b/config_src/drivers/solo_driver/MESO_surface_forcing.F90
index 300c736802..9c67be4829 100644
--- a/config_src/drivers/solo_driver/MESO_surface_forcing.F90
+++ b/config_src/drivers/solo_driver/MESO_surface_forcing.F90
@@ -172,7 +172,7 @@ subroutine MESO_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
       do j=js,je ; do i=is,ie
         !   Set Temp_restore and Salin_restore to the temperature (in degC) and
         ! salinity (in ppt or PSU) that are being restored toward.
-        if (G%mask2dT(i,j) > 0) then
+        if (G%mask2dT(i,j) > 0.0) then
           fluxes%heat_added(i,j) = G%mask2dT(i,j) * &
               ((CS%T_Restore(i,j) - sfc_state%SST(i,j)) * rhoXcp * CS%Flux_const)
           fluxes%vprec(i,j) = - (CS%Rho0 * CS%Flux_const) * &
diff --git a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 b/config_src/drivers/solo_driver/MOM_surface_forcing.F90
index 7c26c2f194..7bca03ca5d 100644
--- a/config_src/drivers/solo_driver/MOM_surface_forcing.F90
+++ b/config_src/drivers/solo_driver/MOM_surface_forcing.F90
@@ -1125,7 +1125,7 @@ subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, US, CS)
 
     if (CS%use_temperature) then
       do j=js,je ; do i=is,ie
-        if (G%mask2dT(i,j) > 0) then
+        if (G%mask2dT(i,j) > 0.0) then
           fluxes%heat_added(i,j) = G%mask2dT(i,j) * &
               ((CS%T_Restore(i,j) - sfc_state%SST(i,j)) * rhoXcp * CS%Flux_const_T)
           fluxes%vprec(i,j) = - (CS%Rho0*CS%Flux_const_S) * &
@@ -1138,7 +1138,7 @@ subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, US, CS)
       enddo ; enddo
     else
       do j=js,je ; do i=is,ie
-        if (G%mask2dT(i,j) > 0) then
+        if (G%mask2dT(i,j) > 0.0) then
           fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * &
                              (CS%G_Earth * CS%Flux_const / CS%Rho0)
         else
@@ -1231,7 +1231,7 @@ subroutine buoyancy_forcing_from_data_override(sfc_state, fluxes, day, dt, G, US
   if (CS%restorebuoy) then
     if (CS%use_temperature) then
       do j=js,je ; do i=is,ie
-        if (G%mask2dT(i,j) > 0) then
+        if (G%mask2dT(i,j) > 0.0) then
           fluxes%heat_added(i,j) = G%mask2dT(i,j) * &
               ((CS%T_Restore(i,j) - sfc_state%SST(i,j)) * rhoXcp * CS%Flux_const_T)
           fluxes%vprec(i,j) = - (CS%Rho0*CS%Flux_const_S) * &
@@ -1244,7 +1244,7 @@ subroutine buoyancy_forcing_from_data_override(sfc_state, fluxes, day, dt, G, US
       enddo ; enddo
     else
       do j=js,je ; do i=is,ie
-        if (G%mask2dT(i,j) > 0) then
+        if (G%mask2dT(i,j) > 0.0) then
           fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * &
                              (CS%G_Earth * CS%Flux_const / CS%Rho0)
         else
@@ -1431,7 +1431,7 @@ subroutine buoyancy_forcing_linear(sfc_state, fluxes, day, dt, G, US, CS)
         y = (G%geoLatCu(I,j)-CS%South_lat)/CS%len_lat
         T_restore = CS%T_south + (CS%T_north-CS%T_south)*y
         S_restore = CS%S_south + (CS%S_north-CS%S_south)*y
-        if (G%mask2dT(i,j) > 0) then
+        if (G%mask2dT(i,j) > 0.0) then
           fluxes%heat_added(i,j) = G%mask2dT(i,j) * &
               ((T_Restore - sfc_state%SST(i,j)) * ((CS%Rho0 * fluxes%C_p) * CS%Flux_const))
           fluxes%vprec(i,j) = - (CS%Rho0*CS%Flux_const) * &
@@ -1446,7 +1446,7 @@ subroutine buoyancy_forcing_linear(sfc_state, fluxes, day, dt, G, US, CS)
       call MOM_error(FATAL, "buoyancy_forcing_linear in MOM_surface_forcing: "// &
                      "RESTOREBUOY to linear not written yet.")
      !do j=js,je ; do i=is,ie
-     !  if (G%mask2dT(i,j) > 0) then
+     !  if (G%mask2dT(i,j) > 0.0) then
      !    fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * &
      !                       (CS%G_Earth * CS%Flux_const / CS%Rho0)
      !  else
diff --git a/src/ALE/coord_adapt.F90 b/src/ALE/coord_adapt.F90
index fe3864fc7a..e5b33103ef 100644
--- a/src/ALE/coord_adapt.F90
+++ b/src/ALE/coord_adapt.F90
@@ -154,7 +154,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex
   ! TODO: this needs to be adjusted to account for vanished layers near topography
 
   ! up (j-1)
-  if (G%mask2dT(i,j-1) > 0.) then
+  if (G%mask2dT(i,j-1) > 0.0) then
     call calculate_density_derivs( &
          0.5 * (tInt(i,j,2:nz) + tInt(i,j-1,2:nz)), &
          0.5 * (sInt(i,j,2:nz) + sInt(i,j-1,2:nz)), &
@@ -166,7 +166,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex
           beta(2:nz)  * (sInt(i,j-1,2:nz) - sInt(i,j,2:nz)))
   endif
   ! down (j+1)
-  if (G%mask2dT(i,j+1) > 0.) then
+  if (G%mask2dT(i,j+1) > 0.0) then
     call calculate_density_derivs( &
          0.5 * (tInt(i,j,2:nz) + tInt(i,j+1,2:nz)), &
          0.5 * (sInt(i,j,2:nz) + sInt(i,j+1,2:nz)), &
@@ -178,7 +178,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex
           beta(2:nz)  * (sInt(i,j+1,2:nz) - sInt(i,j,2:nz)))
   endif
   ! left (i-1)
-  if (G%mask2dT(i-1,j) > 0.) then
+  if (G%mask2dT(i-1,j) > 0.0) then
     call calculate_density_derivs( &
          0.5 * (tInt(i,j,2:nz) + tInt(i-1,j,2:nz)), &
          0.5 * (sInt(i,j,2:nz) + sInt(i-1,j,2:nz)), &
@@ -190,7 +190,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex
           beta(2:nz)  * (sInt(i-1,j,2:nz) - sInt(i,j,2:nz)))
   endif
   ! right (i+1)
-  if (G%mask2dT(i+1,j) > 0.) then
+  if (G%mask2dT(i+1,j) > 0.0) then
     call calculate_density_derivs( &
          0.5 * (tInt(i,j,2:nz) + tInt(i+1,j,2:nz)), &
          0.5 * (sInt(i,j,2:nz) + sInt(i+1,j,2:nz)), &
diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
index ccd70cd9da..6de46f7738 100644
--- a/src/core/MOM_forcing_type.F90
+++ b/src/core/MOM_forcing_type.F90
@@ -2184,12 +2184,12 @@ subroutine set_derived_forcing_fields(forces, fluxes, G, US, Rho0)
       associated(fluxes%ustar_gustless)) then
     do j=js,je ; do i=is,ie
       taux2 = 0.0
-      if ((G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) > 0) &
+      if ((G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) > 0.0) &
         taux2 = (G%mask2dCu(I-1,j) * forces%taux(I-1,j)**2 + &
                  G%mask2dCu(I,j) * forces%taux(I,j)**2) / &
                 (G%mask2dCu(I-1,j) + G%mask2dCu(I,j))
       tauy2 = 0.0
-      if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0) &
+      if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0.0) &
         tauy2 = (G%mask2dCv(i,J-1) * forces%tauy(i,J-1)**2 + &
                  G%mask2dCv(i,J) * forces%tauy(i,J)**2) / &
                 (G%mask2dCv(i,J-1) + G%mask2dCv(i,J))
@@ -3518,15 +3518,15 @@ subroutine homogenize_mech_forcing(forces, G, US, Rho0, UpdateUstar)
   if (do_stress) then
     tx_mean = global_area_mean_u(forces%taux, G, tmp_scale=US%Z_to_L*US%RL2_T2_to_Pa)
     do j=js,je ; do i=isB,ieB
-      if (G%mask2dCu(I,j) > 0.) forces%taux(I,j) = tx_mean
+      if (G%mask2dCu(I,j) > 0.0) forces%taux(I,j) = tx_mean
     enddo ; enddo
     ty_mean = global_area_mean_v(forces%tauy, G, tmp_scale=US%Z_to_L*US%RL2_T2_to_Pa)
     do j=jsB,jeB ; do i=is,ie
-      if (G%mask2dCv(i,J) > 0.) forces%tauy(i,J) = ty_mean
+      if (G%mask2dCv(i,J) > 0.0) forces%tauy(i,J) = ty_mean
     enddo ; enddo
     if (tau2ustar) then
       do j=js,je ; do i=is,ie
-        if (G%mask2dT(i,j) > 0.) forces%ustar(i,j) = sqrt(sqrt(tx_mean**2 + ty_mean**2)*Irho0)
+        if (G%mask2dT(i,j) > 0.0) forces%ustar(i,j) = sqrt(sqrt(tx_mean**2 + ty_mean**2)*Irho0)
       enddo ; enddo
     else
       call homogenize_field_t(forces%ustar, G, tmp_scale=US%Z_to_m*US%s_to_T)
@@ -3689,7 +3689,7 @@ subroutine homogenize_field_t(var, G, tmp_scale)
 
   avg = global_area_mean(var, G, tmp_scale=tmp_scale)
   do j=js,je ; do i=is,ie
-    if (G%mask2dT(i,j) > 0.) var(i,j) = avg
+    if (G%mask2dT(i,j) > 0.0) var(i,j) = avg
   enddo ; enddo
 
 end subroutine homogenize_field_t
@@ -3706,7 +3706,7 @@ subroutine homogenize_field_v(var, G, tmp_scale)
 
   avg = global_area_mean_v(var, G, tmp_scale=tmp_scale)
   do J=jsB,jeB ; do i=is,ie
-    if (G%mask2dCv(i,J) > 0.) var(i,J) = avg
+    if (G%mask2dCv(i,J) > 0.0) var(i,J) = avg
   enddo ; enddo
 
 end subroutine homogenize_field_v
@@ -3723,7 +3723,7 @@ subroutine homogenize_field_u(var, G, tmp_scale)
 
   avg = global_area_mean_u(var, G, tmp_scale=tmp_scale)
   do j=js,je ; do I=isB,ieB
-    if (G%mask2dCu(I,j) > 0.) var(I,j) = avg
+    if (G%mask2dCu(I,j) > 0.0) var(I,j) = avg
   enddo ; enddo
 
 end subroutine homogenize_field_u
diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90
index 0f3b19f2e1..50933abe07 100644
--- a/src/diagnostics/MOM_wave_speed.F90
+++ b/src/diagnostics/MOM_wave_speed.F90
@@ -262,7 +262,7 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_
     endif
 
     ! From this point, we can work on individual columns without causing memory to have page faults.
-    do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+    do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
       if (use_EOS) then
         pres(1) = 0.0 ; H_top(1) = 0.0
         do K=2,kf(i)
@@ -815,7 +815,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos)
 
     ! From this point, we can work on individual columns without causing memory to have page faults.
     do i=is,ie
-      if (G%mask2dT(i,j) > 0.5) then
+      if (G%mask2dT(i,j) > 0.0) then
         if (use_EOS) then
           pres(1) = 0.0 ; H_top(1) = 0.0
           do K=2,kf(i)
diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90
index c833e973c5..81c8799828 100644
--- a/src/diagnostics/MOM_wave_structure.F90
+++ b/src/diagnostics/MOM_wave_structure.F90
@@ -272,12 +272,12 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
 
     ! From this point, we can work on individual columns without causing memory
     ! to have page faults.
-    do i=is,ie ; if (cn(i,j)>0.0)then
+    do i=is,ie ; if (cn(i,j) > 0.0) then
       !----for debugging, remove later----
       ig = i + G%idg_offset ; jg = j + G%jdg_offset
       !if (ig == CS%int_tide_source_x .and. jg == CS%int_tide_source_y) then
       !-----------------------------------
-      if (G%mask2dT(i,j) > 0.5) then
+      if (G%mask2dT(i,j) > 0.0) then
 
         gprime(:) = 0.0 ! init gprime
         pres(:) = 0.0 ! init pres
@@ -567,7 +567,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
       !else     ! if at test point - delete later
       !  return ! if at test point - delete later
       !endif    ! if at test point - delete later
-      endif ! mask2dT > 0.5?
+      endif ! mask2dT > 0.0?
     else
       ! if cn=0.0, default to zero
       nzm                       = nz+1! could use actual values
diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90
index f711ef3f87..446cdb8b45 100644
--- a/src/parameterizations/vertical/MOM_ALE_sponge.F90
+++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90
@@ -230,7 +230,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h,
   ! number of columns to be restored
   CS%num_col = 0 ; CS%fldno = 0
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
-    if ((Iresttime(i,j)>0.0) .and. (G%mask2dT(i,j)>0)) &
+    if ((Iresttime(i,j) > 0.0) .and. (G%mask2dT(i,j) > 0.0)) &
       CS%num_col = CS%num_col + 1
   enddo ; enddo
 
@@ -241,7 +241,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h,
     ! pass indices, restoring time to the CS structure
     col = 1
     do j=G%jsc,G%jec ; do i=G%isc,G%iec
-      if ((Iresttime(i,j)>0.0) .and. (G%mask2dT(i,j)>0)) then
+      if ((Iresttime(i,j) > 0.0) .and. (G%mask2dT(i,j) > 0.0)) then
         CS%col_i(col) = i ; CS%col_j(col) = j
         CS%Iresttime_col(col) = Iresttime(i,j)
         col = col +1
@@ -282,7 +282,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h,
       enddo ; enddo
     endif
     do j=G%jsc,G%jec ; do I=G%iscB,G%iecB
-       if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) &
+       if ((Iresttime_u(I,j) > 0.0) .and. (G%mask2dCu(I,j) > 0.0)) &
           CS%num_col_u = CS%num_col_u + 1
     enddo ; enddo
 
@@ -295,7 +295,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h,
       ! Store the column indices and restoring rates in the CS structure
       col = 1
       do j=G%jsc,G%jec ; do I=G%iscB,G%iecB
-        if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then
+        if ((Iresttime_u(I,j) > 0.0) .and. (G%mask2dCu(I,j) > 0.0)) then
           CS%col_i_u(col) = I ; CS%col_j_u(col) = j
           CS%Iresttime_col_u(col) = Iresttime_u(I,j)
           col = col + 1
@@ -326,7 +326,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h,
       enddo ; enddo
     endif
     do J=G%jscB,G%jecB; do i=G%isc,G%iec
-      if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) &
+      if ((Iresttime_v(i,J) > 0.0) .and. (G%mask2dCv(i,J) > 0.0)) &
         CS%num_col_v = CS%num_col_v + 1
     enddo ; enddo
 
@@ -339,7 +339,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h,
       ! pass indices, restoring time to the CS structure
       col = 1
       do J=G%jscB,G%jecB ; do i=G%isc,G%iec
-        if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) then
+        if ((Iresttime_v(i,J) > 0.0) .and. (G%mask2dCv(i,J) > 0.0)) then
           CS%col_i_v(col) = i ; CS%col_j_v(col) = j
           CS%Iresttime_col_v(col) = Iresttime_v(i,j)
           col = col + 1
@@ -492,7 +492,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest
   ! number of columns to be restored
   CS%num_col = 0 ; CS%fldno = 0
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
-    if ((Iresttime(i,j)>0.0) .and. (G%mask2dT(i,j)>0)) &
+    if ((Iresttime(i,j) > 0.0) .and. (G%mask2dT(i,j) > 0.0)) &
       CS%num_col = CS%num_col + 1
   enddo ; enddo
   if (CS%num_col > 0) then
@@ -502,7 +502,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest
     ! pass indices, restoring time to the CS structure
     col = 1
     do j=G%jsc,G%jec ; do i=G%isc,G%iec
-      if ((Iresttime(i,j)>0.0) .and. (G%mask2dT(i,j)>0)) then
+      if ((Iresttime(i,j) > 0.0) .and. (G%mask2dT(i,j) > 0.0)) then
         CS%col_i(col) = i ; CS%col_j(col) = j
         CS%Iresttime_col(col) = Iresttime(i,j)
         col = col + 1
@@ -532,7 +532,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest
     endif
     CS%num_col_u = 0 ;
     do j=G%jsc,G%jec; do I=G%iscB,G%iecB
-      if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) &
+      if ((Iresttime_u(I,j) > 0.0) .and. (G%mask2dCu(I,j) > 0.0)) &
         CS%num_col_u = CS%num_col_u + 1
     enddo ; enddo
     if (CS%num_col_u > 0) then
@@ -542,7 +542,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest
       ! pass indices, restoring time to the CS structure
       col = 1
       do j=G%jsc,G%jec ; do I=G%iscB,G%iecB
-        if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then
+        if ((Iresttime_u(I,j) > 0.0) .and. (G%mask2dCu(I,j) > 0.0)) then
           CS%col_i_u(col) = i ; CS%col_j_u(col) = j
           CS%Iresttime_col_u(col) = Iresttime_u(i,j)
           col = col + 1
@@ -564,7 +564,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest
     endif
     CS%num_col_v = 0 ;
     do J=G%jscB,G%jecB; do i=G%isc,G%iec
-      if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) &
+      if ((Iresttime_v(i,J) > 0.0) .and. (G%mask2dCv(i,J) > 0.0)) &
         CS%num_col_v = CS%num_col_v + 1
     enddo ; enddo
     if (CS%num_col_v > 0) then
@@ -574,7 +574,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest
       ! pass indices, restoring time to the CS structure
       col = 1
       do J=G%jscB,G%jecB ; do i=G%isc,G%iec
-        if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) then
+        if ((Iresttime_v(i,J) > 0.0) .and. (G%mask2dCv(i,J) > 0.0)) then
           CS%col_i_v(col) = i ; CS%col_j_v(col) = j
           CS%Iresttime_col_v(col) = Iresttime_v(i,j)
           col = col + 1
diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90
index 6f1988c295..7aec4c0331 100644
--- a/src/parameterizations/vertical/MOM_diabatic_aux.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90
@@ -527,7 +527,7 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb, zero_mix)
   do j=js,je
     do i=is,ie
       sum_area = G%areaCu(I-1,j) + G%areaCu(I,j)
-      if (sum_area>0.0) then
+      if (sum_area > 0.0) then
         Idenom = sqrt(0.5*G%IareaT(i,j) / sum_area)
         a_w(i) = G%areaCu(I-1,j) * Idenom
         a_e(i) = G%areaCu(I,j) * Idenom
@@ -536,7 +536,7 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb, zero_mix)
       endif
 
       sum_area = G%areaCv(i,J-1) + G%areaCv(i,J)
-      if (sum_area>0.0) then
+      if (sum_area > 0.0) then
         Idenom = sqrt(0.5*G%IareaT(i,j) / sum_area)
         a_s(i) = G%areaCv(i,J-1) * Idenom
         a_n(i) = G%areaCv(i,J) * Idenom
@@ -618,7 +618,7 @@ subroutine set_pen_shortwave(optics, fluxes, G, GV, US, CS, opacity, tracer_flow
       ! same value is assumed for all layers.
       call time_interp_external(CS%sbc_chl, CS%Time, chl_2d)
       do j=js,je ; do i=is,ie
-        if ((G%mask2dT(i,j) > 0.5) .and. (chl_2d(i,j) < 0.0)) then
+        if ((G%mask2dT(i,j) > 0.0) .and. (chl_2d(i,j) < 0.0)) then
           write(mesg,'(" Time_interp negative chl of ",(1pe12.4)," at i,j = ",&
                     & 2(i3), "lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
                      chl_2d(i,j), i, j, G%geoLonT(i,j), G%geoLatT(i,j)
@@ -768,7 +768,7 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US,
       ! endif
       call calculate_density(T_subML, S_subML, pRef_N2, rho_subML, tv%eqn_of_state, EOSdom)
       call calculate_density(T_deeper, S_deeper, pRef_N2, rho_deeper, tv%eqn_of_state, EOSdom)
-      do i=is,ie ; if ((G%mask2dT(i,j)>0.5) .and. N2_region_set(i)) then
+      do i=is,ie ; if ((G%mask2dT(i,j) > 0.0) .and. N2_region_set(i)) then
         subMLN2(i,j) =  gE_rho0 * (rho_deeper(i) - rho_subML(i)) / (GV%H_to_z * dH_N2(i))
       endif ; enddo
     endif
@@ -1227,7 +1227,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
       else
         netMassIn(i) = netMassInOut(i) - netMassOut(i)
       endif
-      if (G%mask2dT(i,j)>0.0) then
+      if (G%mask2dT(i,j) > 0.0) then
         fluxes%netMassOut(i,j) = netMassOut(i)
         fluxes%netMassIn(i,j) = netMassIn(i)
       else
diff --git a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90
index 379275196e..eb592fb890 100644
--- a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90
+++ b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90
@@ -81,7 +81,7 @@ subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV, US, CS, Kd_int)
          "Module must be initialized before it is used.")
 
 !$OMP do
-  do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+  do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
     if (present(Kd_int) .and. .not.CS%use_test_Kh_profile) then
       do k=1,nz+1 ; Kd(K) = CS%test_Kh_scaling*Kd_int(i,j,K) ; enddo
     else
diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90
index 99dd38135d..98ab99616f 100644
--- a/src/parameterizations/vertical/MOM_energetic_PBL.F90
+++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90
@@ -390,7 +390,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS
     ! homogenizing the shortwave heating within that cell.  This sets the energy
     ! and ustar and wstar available to drive mixing at the first interior
     ! interface.
-    do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+    do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
 
       ! Copy the thicknesses and other fields to 1-d arrays.
       do k=1,nz
diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90
index fdc38ebf1e..45d442f98c 100644
--- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90
+++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90
@@ -288,7 +288,7 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
       enddo ; enddo
     endif
 
-    do i=is,ie ; do_i(i) = (G%mask2dT(i,j) > 0.5) ; enddo
+    do i=is,ie ; do_i(i) = (G%mask2dT(i,j) > 0.0) ; enddo
     do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; enddo
     do i=is,ie ; dsp1_ds(i,nz) = 0.0 ; enddo
 
diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90
index 300cdcbe1e..7af12aee6f 100644
--- a/src/parameterizations/vertical/MOM_internal_tide_input.F90
+++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90
@@ -253,7 +253,7 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot)
     do i=is,ie
       hb(i) = 0.0 ; dRho_bot(i) = 0.0
       z_from_bot(i) = 0.5*GV%H_to_Z*h(i,j,nz)
-      do_i(i) = (G%mask2dT(i,j) > 0.5)
+      do_i(i) = (G%mask2dT(i,j) > 0.0)
       h_amp(i) = sqrt(h2(i,j))
     enddo
 
diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90
index a44a7aee95..8af8727246 100644
--- a/src/parameterizations/vertical/MOM_kappa_shear.F90
+++ b/src/parameterizations/vertical/MOM_kappa_shear.F90
@@ -197,7 +197,7 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, &
 !---------------------------------------
 ! Work on each column.
 !---------------------------------------
-    do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+    do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
     ! call cpu_clock_begin(id_clock_setup)
       ! Store a transposed version of the initial arrays.
       ! Any elimination of massless layers would occur here.
diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90
index 9aa8fafd14..66b7a7746e 100644
--- a/src/parameterizations/vertical/MOM_opacity.F90
+++ b/src/parameterizations/vertical/MOM_opacity.F90
@@ -290,7 +290,7 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir
   if (present(chl_3d)) then
     do j=js,je ; do i=is,ie ; chl_data(i,j) = chl_3d(i,j,1) ; enddo ; enddo
     do k=1,nz ; do j=js,je ; do i=is,ie
-      if ((G%mask2dT(i,j) > 0.5) .and. (chl_3d(i,j,k) < 0.0)) then
+      if ((G%mask2dT(i,j) > 0.0) .and. (chl_3d(i,j,k) < 0.0)) then
         write(mesg,'(" Negative chl_3d of ",(1pe12.4)," found at i,j,k = ", &
                   & 3(1x,i3), " lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
                    chl_3d(i,j,k), i, j, k, G%geoLonT(i,j), G%geoLatT(i,j)
@@ -300,7 +300,7 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir
   elseif (present(chl_2d)) then
     do j=js,je ; do i=is,ie ; chl_data(i,j) = chl_2d(i,j) ; enddo ; enddo
     do j=js,je ; do i=is,ie
-      if ((G%mask2dT(i,j) > 0.5) .and. (chl_2d(i,j) < 0.0)) then
+      if ((G%mask2dT(i,j) > 0.0) .and. (chl_2d(i,j) < 0.0)) then
         write(mesg,'(" Negative chl_2d of ",(1pe12.4)," at i,j = ", &
                   & 2(i3), "lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
                    chl_data(i,j), i, j, G%geoLonT(i,j), G%geoLatT(i,j)
@@ -316,7 +316,7 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir
       !$OMP parallel do default(shared) private(SW_vis_tot,SW_nir_tot)
       do j=js,je ; do i=is,ie
         SW_vis_tot = 0.0 ; SW_nir_tot = 0.0
-        if (G%mask2dT(i,j) > 0.5) then
+        if (G%mask2dT(i,j) > 0.0) then
           if (multiband_vis_input) then
             SW_vis_tot = sw_vis_dir(i,j) + sw_vis_dif(i,j)
           elseif (total_sw_input) then
@@ -344,7 +344,7 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir
       !$OMP parallel do default(shared) private(SW_pen_tot)
       do j=js,je ; do i=is,ie
         SW_pen_tot = 0.0
-        if (G%mask2dT(i,j) > 0.5) then
+        if (G%mask2dT(i,j) > 0.0) then
           if (multiband_vis_input) then
             SW_pen_tot = SW_pen_frac_morel(chl_data(i,j)) * (sw_vis_dir(i,j) + sw_vis_dif(i,j))
           elseif (total_sw_input) then
@@ -385,7 +385,7 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir
       case (MOREL_88)
         do j=js,je ; do i=is,ie
           optics%opacity_band(1,i,j,k) = CS%opacity_land_value
-          if (G%mask2dT(i,j) > 0.5) &
+          if (G%mask2dT(i,j) > 0.0) &
             optics%opacity_band(1,i,j,k) = US%Z_to_m * opacity_morel(chl_data(i,j))
 
           do n=2,optics%nbands
diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90
index 061fe776e1..0d5f1a9ca9 100644
--- a/src/parameterizations/vertical/MOM_set_diffusivity.F90
+++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90
@@ -503,7 +503,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
 
 
     ! This adds the diffusion sustained by the energy extracted from the flow by the bottom drag.
-    if (CS%bottomdraglaw .and. (CS%BBL_effic>0.0)) then
+    if (CS%bottomdraglaw .and. (CS%BBL_effic > 0.0)) then
       if (CS%use_LOTW_BBL_diffusivity) then
         call add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int_2d, G, GV, US, CS, &
                                       dd%Kd_BBL, Kd_lay_2d)
@@ -812,7 +812,7 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, &
 
   do i=is,ie
     htot(i) = GV%H_to_Z*(h(i,j,nz) - GV%Angstrom_H) ; maxEnt(i,nz) = 0.0
-    do_i(i) = (G%mask2dT(i,j) > 0.5)
+    do_i(i) = (G%mask2dT(i,j) > 0.0)
   enddo
   do k=nz-1,kb_min,-1
     i_rem = 0
@@ -969,7 +969,7 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, &
   do i=is,ie
     hb(i) = 0.0 ; dRho_bot(i) = 0.0 ; h_amp(i) = 0.0
     z_from_bot(i) = 0.5*GV%H_to_Z*h(i,j,nz)
-    do_i(i) = (G%mask2dT(i,j) > 0.5)
+    do_i(i) = (G%mask2dT(i,j) > 0.0)
   enddo
   if (CS%use_tidal_mixing) call tidal_mixing_h_amp(h_amp, G, j, CS%tidal_mixing)
 
@@ -1001,7 +1001,7 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, &
       N2_bot(i) = (G_Rho0 * drho_bot(i)) / hb(i)
     else ;  N2_bot(i) = 0.0 ; endif
     z_from_bot(i) = 0.5*GV%H_to_Z*h(i,j,nz)
-    do_i(i) = (G%mask2dT(i,j) > 0.5)
+    do_i(i) = (G%mask2dT(i,j) > 0.0)
   enddo
 
   do k=nz,2,-1
@@ -1195,7 +1195,7 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, &
 
   do_diag_Kd_BBL = associated(Kd_BBL)
 
-  if (.not.(CS%bottomdraglaw .and. (CS%BBL_effic>0.0))) return
+  if (.not.(CS%bottomdraglaw .and. (CS%BBL_effic > 0.0))) return
 
   cdrag_sqrt = sqrt(CS%cdrag)
   TKE_Ray = 0.0 ; Rayleigh_drag = .false.
@@ -1238,7 +1238,7 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, &
 
     gh_sum_top(i) = R0_g * 400.0 * ustar_h**2
 
-    do_i(i) = (G%mask2dT(i,j) > 0.5)
+    do_i(i) = (G%mask2dT(i,j) > 0.0)
     htot(i) = GV%H_to_Z*h(i,j,nz)
     rho_htot(i) = GV%Rlay(nz)*(GV%H_to_Z*h(i,j,nz))
     Rho_top(i) = GV%Rlay(1)
@@ -1259,7 +1259,7 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, &
     if (.not.domore) exit
   enddo ! k-loop
 
-  do i=is,ie ; do_i(i) = (G%mask2dT(i,j) > 0.5) ; enddo
+  do i=is,ie ; do_i(i) = (G%mask2dT(i,j) > 0.0) ; enddo
   do k=nz-1,kb_min,-1
     i_rem = 0
     do i=is,ie ; if (do_i(i)) then
@@ -1409,7 +1409,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int
   real, parameter :: von_karm = 0.41 ! Von Karman constant (http://en.wikipedia.org/wiki/Von_Karman_constant)
   logical :: do_diag_Kd_BBL
 
-  if (.not.(CS%bottomdraglaw .and. (CS%BBL_effic>0.0))) return
+  if (.not.(CS%bottomdraglaw .and. (CS%BBL_effic > 0.0))) return
   do_diag_Kd_BBL = associated(Kd_BBL)
 
   N2_min = 0.
@@ -1568,7 +1568,7 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, Kd_int, G, GV, US, CS, TKE_to_Kd,
 
   if (.not.CS%ML_radiation) return
 
-  do i=is,ie ; h_ml(i) = 0.0 ; do_i(i) = (G%mask2dT(i,j) > 0.5) ; enddo
+  do i=is,ie ; h_ml(i) = 0.0 ; do_i(i) = (G%mask2dT(i,j) > 0.0) ; enddo
   do k=1,kml ; do i=is,ie ; h_ml(i) = h_ml(i) + GV%H_to_Z*h(i,j,k) ; enddo ; enddo
 
   do i=is,ie ; if (do_i(i)) then
@@ -1734,7 +1734,7 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC)
     ! Determine ustar and the square magnitude of the velocity in the
     ! bottom boundary layer. Together these give the TKE source and
     ! vertical decay scale.
-    do i=is,ie ; if ((G%mask2dCv(i,J) > 0.5) .and. (cdrag_sqrt*visc%bbl_thick_v(i,J) > 0.0)) then
+    do i=is,ie ; if ((G%mask2dCv(i,J) > 0.0) .and. (cdrag_sqrt*visc%bbl_thick_v(i,J) > 0.0)) then
       do_i(i) = .true. ; vhtot(i) = 0.0 ; htot(i) = 0.0
       vstar(i,J) = visc%Kv_bbl_v(i,J) / (cdrag_sqrt*visc%bbl_thick_v(i,J))
     else
@@ -1775,7 +1775,7 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC)
       endif ; enddo
       if (.not.domore) exit
     enddo
-    do i=is,ie ; if ((G%mask2dCv(i,J) > 0.5) .and. (htot(i) > 0.0)) then
+    do i=is,ie ; if ((G%mask2dCv(i,J) > 0.0) .and. (htot(i) > 0.0)) then
       v2_bbl(i,J) = (vhtot(i)*vhtot(i))/(htot(i)*htot(i))
     else
       v2_bbl(i,J) = 0.0
@@ -1783,7 +1783,7 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC)
   enddo
   !$OMP do
   do j=js,je
-    do I=is-1,ie ; if ((G%mask2dCu(I,j) > 0.5) .and. (cdrag_sqrt*visc%bbl_thick_u(I,j) > 0.0))  then
+    do I=is-1,ie ; if ((G%mask2dCu(I,j) > 0.0) .and. (cdrag_sqrt*visc%bbl_thick_u(I,j) > 0.0))  then
       do_i(I) = .true. ; uhtot(I) = 0.0 ; htot(I) = 0.0
       ustar(I) = visc%Kv_bbl_u(I,j) / (cdrag_sqrt*visc%bbl_thick_u(I,j))
     else
@@ -1823,7 +1823,7 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC)
       endif ; enddo
       if (.not.domore) exit
     enddo
-    do I=is-1,ie ; if ((G%mask2dCu(I,j) > 0.5) .and. (htot(i) > 0.0)) then
+    do I=is-1,ie ; if ((G%mask2dCu(I,j) > 0.0) .and. (htot(i) > 0.0)) then
       u2_bbl(I) = (uhtot(I)*uhtot(I))/(htot(I)*htot(I))
     else
       u2_bbl(I) = 0.0
diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90
index 92e7b44128..93a71f6b64 100644
--- a/src/parameterizations/vertical/MOM_set_viscosity.F90
+++ b/src/parameterizations/vertical/MOM_set_viscosity.F90
@@ -385,15 +385,13 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
       if (j<G%Jsc) cycle
       is = Isq ; ie = Ieq
       do i=is,ie
-        do_i(i) = .false.
-        if (G%mask2dCu(I,j) > 0) do_i(i) = .true.
+        do_i(i) = (G%mask2dCu(I,j) > 0.0)
       enddo
     else
       ! m=2 refers to v-points
       is = G%isc ; ie = G%iec
       do i=is,ie
-        do_i(i) = .false.
-        if (G%mask2dCv(i,J) > 0) do_i(i) = .true.
+        do_i(i) = (G%mask2dCv(i,J) > 0.0)
       enddo
     endif
 
diff --git a/src/parameterizations/vertical/MOM_sponge.F90 b/src/parameterizations/vertical/MOM_sponge.F90
index d0d64079c3..dee80b9e40 100644
--- a/src/parameterizations/vertical/MOM_sponge.F90
+++ b/src/parameterizations/vertical/MOM_sponge.F90
@@ -131,7 +131,7 @@ subroutine initialize_sponge(Iresttime, int_height, G, param_file, CS, GV, &
 
   CS%num_col = 0 ; CS%fldno = 0
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
-    if ((Iresttime(i,j)>0.0) .and. (G%mask2dT(i,j)>0)) &
+    if ((Iresttime(i,j) > 0.0) .and. (G%mask2dT(i,j) > 0.0)) &
       CS%num_col = CS%num_col + 1
   enddo ; enddo
 
@@ -143,7 +143,7 @@ subroutine initialize_sponge(Iresttime, int_height, G, param_file, CS, GV, &
 
     col = 1
     do j=G%jsc,G%jec ; do i=G%isc,G%iec
-      if ((Iresttime(i,j)>0.0) .and. (G%mask2dT(i,j)>0)) then
+      if ((Iresttime(i,j) > 0.0) .and. (G%mask2dT(i,j) > 0.0)) then
         CS%col_i(col) = i ; CS%col_j(col) = j
         CS%Iresttime_col(col) = Iresttime(i,j)
         col = col +1
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index d384500c3d..d6c93a5211 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -256,7 +256,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
   !$OMP                 private(do_i,surface_stress,zDS,stress,h_a,hfr, &
   !$OMP                         b_denom_1,b1,d1,c1)
   do j=G%jsc,G%jec
-    do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0) ; enddo
+    do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0.0) ; enddo
 
     ! When mixing down Eulerian current + Stokes drift add before calling solver
     if (DoStokesMixing) then ; do k=1,nz ; do I=Isq,Ieq
@@ -383,7 +383,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
   !$OMP               private(do_i,surface_stress,zDS,stress,h_a,hfr, &
   !$OMP                       b_denom_1,b1,d1,c1)
   do J=Jsq,Jeq
-    do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0) ; enddo
+    do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0.0) ; enddo
 
     ! When mixing down Eulerian current + Stokes drift add before calling solver
     if (DoStokesMixing) then ; do k=1,nz ; do i=is,ie
@@ -593,7 +593,7 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS)
   ! Find the zonal viscous remnant using a modification of a standard tridagonal solver.
   !$OMP parallel do default(shared) firstprivate(Ray) private(do_i,b_denom_1,b1,d1,c1)
   do j=G%jsc,G%jec
-    do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0) ; enddo
+    do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0.0) ; enddo
 
     if (CS%Channel_drag) then ; do k=1,nz ; do I=Isq,Ieq
       Ray(I,k) = visc%Ray_u(I,j,k)
@@ -622,7 +622,7 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS)
   ! Now find the meridional viscous remnant using the robust tridiagonal solver.
   !$OMP parallel do default(shared) firstprivate(Ray) private(do_i,b_denom_1,b1,d1,c1)
   do J=Jsq,Jeq
-    do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0) ; enddo
+    do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0.0) ; enddo
 
     if (CS%Channel_drag) then ; do k=1,nz ; do i=is,ie
       Ray(i,k) = visc%Ray_v(i,J,k)
@@ -766,7 +766,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
   !$OMP                                     OBC,h_neglect,dt,I_valBL,Kv_u,a_cpl_max) &
   !$OMP                     firstprivate(i_hbbl)
   do j=G%Jsc,G%Jec
-    do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0) ; enddo
+    do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0.0) ; enddo
 
     if (CS%bottomdraglaw) then ; do I=Isq,Ieq
       kv_bbl(I) = visc%Kv_bbl_u(I,j)
@@ -931,7 +931,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
   !$OMP                                  OBC,h_neglect,dt,I_valBL,Kv_v,a_cpl_max) &
   !$OMP                     firstprivate(i_hbbl)
   do J=Jsq,Jeq
-    do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0) ; enddo
+    do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0.0) ; enddo
 
     if (CS%bottomdraglaw) then ; do i=is,ie
       kv_bbl(i) = visc%Kv_bbl_v(i,J)
diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90
index 31acb51160..5e8632f1ca 100644
--- a/src/tracer/MOM_generic_tracer.F90
+++ b/src/tracer/MOM_generic_tracer.F90
@@ -361,7 +361,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file,
     grid_tmask(:,:,:) = 0.0
     grid_kmt(:,:) = 0
     do j = G%jsd, G%jed ; do i = G%isd, G%ied
-      if (G%mask2dT(i,j) > 0) then
+      if (G%mask2dT(i,j) > 0.0) then
         grid_tmask(i,j,:) = 1.0
         grid_kmt(i,j) = GV%ke ! Tell the code that a layer thicker than 1m is the bottom layer.
       endif
diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90
index 61715e044b..08361ba9af 100644
--- a/src/tracer/MOM_neutral_diffusion.F90
+++ b/src/tracer/MOM_neutral_diffusion.F90
@@ -325,7 +325,7 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf)
     call pass_var(hbl,G%Domain)
     ! get k-indices and zeta
     do j=G%jsc-1, G%jec+1 ; do i=G%isc-1,G%iec+1
-      if (G%mask2dT(i,j) > 0.) then
+      if (G%mask2dT(i,j) > 0.0) then
         call boundary_k_range(SURFACE, G%ke, h(i,j,:), hbl(i,j), k_top(i,j), zeta_top(i,j), k_bot(i,j), zeta_bot(i,j))
       endif
     enddo; enddo
@@ -463,7 +463,7 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf)
 
   ! Neutral surface factors at U points
   do j = G%jsc, G%jec ; do I = G%isc-1, G%iec
-    if (G%mask2dCu(I,j) > 0.) then
+    if (G%mask2dCu(I,j) > 0.0) then
       if (CS%continuous_reconstruction) then
         call find_neutral_surface_positions_continuous(GV%ke,                                    &
                 CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:),            &
@@ -484,7 +484,7 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf)
 
   ! Neutral surface factors at V points
   do J = G%jsc-1, G%jec ; do i = G%isc, G%iec
-    if (G%mask2dCv(i,J) > 0.) then
+    if (G%mask2dCv(i,J) > 0.0) then
       if (CS%continuous_reconstruction) then
         call find_neutral_surface_positions_continuous(GV%ke,                                              &
                 CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:),           &
@@ -519,10 +519,10 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf)
       enddo ; enddo ; enddo
     else
       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
+        if (G%mask2dCu(I,j) > 0.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
+        if (G%mask2dCv(i,J) > 0.0) CS%vhEff(i,J,k) = CS%vhEff(i,J,k) * pa_to_H
       enddo ; enddo ; enddo
     endif
   endif
diff --git a/src/tracer/MOM_offline_aux.F90 b/src/tracer/MOM_offline_aux.F90
index bdd6be4fe0..fdf199d71e 100644
--- a/src/tracer/MOM_offline_aux.F90
+++ b/src/tracer/MOM_offline_aux.F90
@@ -179,7 +179,7 @@ subroutine limit_mass_flux_3d(G, GV, uh, vh, ea, eb, h_pre)
                  (max(0.0, -vh(i,J-1,k)) + max(0.0, vh(i,J,k)))) + &
                 (max(0.0, top_flux(i,j,k)) + max(0.0, bottom_flux(i,j,k))) * G%areaT(i,j)
 
-    if (pos_flux>hvol .and. pos_flux>0.0) then
+    if ((pos_flux > hvol) .and. (pos_flux > 0.0)) then
       scale_factor = (hvol / pos_flux) * max_off_cfl
     else ! Don't scale
       scale_factor = 1.0
@@ -405,7 +405,7 @@ subroutine distribute_residual_uh_upwards(G, GV, hvol, uh)
       do k=1,nz
         uh_remain = uh2d(I,k)
         uh2d(I,k) = 0.0
-        if (abs(uh_remain)>0.0) then
+        if (abs(uh_remain) > 0.0) then
           do k_rev = k,1,-1
             uh_sum = uh_remain + uh2d(I,k_rev)
             if (uh_sum<0.0) then ! Transport to the left
@@ -436,7 +436,7 @@ subroutine distribute_residual_uh_upwards(G, GV, hvol, uh)
           enddo ! k_rev
         endif
 
-        if (abs(uh_remain)>0.0) then
+        if (abs(uh_remain) > 0.0) then
           if (k<nz) then
             uh2d(I,k+1) = uh2d(I,k+1) + uh_remain
           else
@@ -503,7 +503,7 @@ subroutine distribute_residual_vh_upwards(G, GV, hvol, vh)
       do k=1,nz
         vh_remain = vh2d(J,k)
         vh2d(J,k) = 0.0
-        if (abs(vh_remain)>0.0) then
+        if (abs(vh_remain) > 0.0) then
           do k_rev = k,1,-1
             vh_sum = vh_remain + vh2d(J,k_rev)
             if (vh_sum<0.0) then ! Transport to the left
@@ -535,7 +535,7 @@ subroutine distribute_residual_vh_upwards(G, GV, hvol, vh)
           enddo ! k_rev
         endif
 
-        if (abs(vh_remain)>0.0) then
+        if (abs(vh_remain) > 0.0) then
          if (k<nz) then
             vh2d(J,k+1) = vh2d(J,k+1) + vh_remain
           else
diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90
index 66adfb3695..a241776d59 100644
--- a/src/tracer/MOM_tracer_diabatic.F90
+++ b/src/tracer/MOM_tracer_diabatic.F90
@@ -147,14 +147,14 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, &
       enddo
 
       ! Now solve the tridiagonal equation for the tracer concentrations.
-      do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         b_denom_1 = h_minus_dsink(i,1) + ea(i,j,1) + h_neglect
         b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
         d1(i) = b_denom_1 * b1(i)
         h_tr = h_old(i,j,1) + h_neglect
         tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j)
       endif ; enddo
-      do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         c1(i,k) = eb(i,j,k-1) * b1(i)
         b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ea(i,j,k) + sink(i,K)) + &
                     h_neglect
@@ -164,7 +164,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, &
         tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + &
                              (ea(i,j,k) + sink(i,K)) * tr(i,j,k-1))
       endif ; enddo ; enddo
-      do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         c1(i,nz) = eb(i,j,nz-1) * b1(i)
         b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ea(i,j,nz) + sink(i,nz)) + &
                     h_neglect
@@ -173,25 +173,25 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, &
         tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + &
                               (ea(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1))
       endif ; enddo
-      if (present(btm_reservoir)) then ; do i=is,ie ; if (G%mask2dT(i,j)>0.5) then
+      if (present(btm_reservoir)) then ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         btm_reservoir(i,j) = btm_reservoir(i,j) + (sink(i,nz+1)*tr(i,j,nz)) * GV%H_to_RZ
       endif ; enddo ; endif
 
-      do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
       endif ; enddo ; enddo
     enddo
   else
     !$OMP do
     do j=js,je
-      do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         h_tr = h_old(i,j,1) + h_neglect
         b_denom_1 = h_tr + ea(i,j,1)
         b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
         d1(i) = h_tr * b1(i)
         tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j)
       endif ; enddo
-      do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         c1(i,k) = eb(i,j,k-1) * b1(i)
         h_tr = h_old(i,j,k) + h_neglect
         b_denom_1 = h_tr + d1(i) * ea(i,j,k)
@@ -199,7 +199,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, &
         d1(i) = b_denom_1 * b1(i)
         tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ea(i,j,k) * tr(i,j,k-1))
       endif ; enddo ; enddo
-      do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         c1(i,nz) = eb(i,j,nz-1) * b1(i)
         h_tr = h_old(i,j,nz) + h_neglect
         b_denom_1 = h_tr + d1(i)*ea(i,j,nz)
@@ -207,7 +207,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, &
         tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + &
                               ea(i,j,nz) * tr(i,j,nz-1))
       endif ; enddo
-      do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
       endif ; enddo ; enddo
     enddo
@@ -345,14 +345,14 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, &
       enddo
 
       ! Now solve the tridiagonal equation for the tracer concentrations.
-      do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         b_denom_1 = h_minus_dsink(i,1) + ent(i,j,1) + h_neglect
         b1(i) = 1.0 / (b_denom_1 + ent(i,j,2))
         d1(i) = b_denom_1 * b1(i)
         h_tr = h_old(i,j,1) + h_neglect
         tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j)
       endif ; enddo
-      do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         c1(i,k) = ent(i,j,K) * b1(i)
         b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ent(i,j,K) + sink(i,K)) + &
                     h_neglect
@@ -362,7 +362,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, &
         tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + &
                              (ent(i,j,K) + sink(i,K)) * tr(i,j,k-1))
       endif ; enddo ; enddo
-      do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         c1(i,nz) = ent(i,j,nz) * b1(i)
         b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ent(i,j,nz) + sink(i,nz)) + &
                     h_neglect
@@ -371,25 +371,25 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, &
         tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + &
                               (ent(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1))
       endif ; enddo
-      if (present(btm_reservoir)) then ; do i=is,ie ; if (G%mask2dT(i,j)>0.5) then
+      if (present(btm_reservoir)) then ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         btm_reservoir(i,j) = btm_reservoir(i,j) + (sink(i,nz+1)*tr(i,j,nz)) * GV%H_to_RZ
       endif ; enddo ; endif
 
-      do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
       endif ; enddo ; enddo
     enddo
   else
     !$OMP do
     do j=js,je
-      do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         h_tr = h_old(i,j,1) + h_neglect
         b_denom_1 = h_tr + ent(i,j,1)
         b1(i) = 1.0 / (b_denom_1 + ent(i,j,2))
         d1(i) = h_tr * b1(i)
         tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j)
       endif ; enddo
-      do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         c1(i,k) = ent(i,j,K) * b1(i)
         h_tr = h_old(i,j,k) + h_neglect
         b_denom_1 = h_tr + d1(i) * ent(i,j,K)
@@ -397,7 +397,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, &
         d1(i) = b_denom_1 * b1(i)
         tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ent(i,j,K) * tr(i,j,k-1))
       endif ; enddo ; enddo
-      do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         c1(i,nz) = ent(i,j,nz) * b1(i)
         h_tr = h_old(i,j,nz) + h_neglect
         b_denom_1 = h_tr + d1(i)*ent(i,j,nz)
@@ -405,7 +405,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, &
         tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + &
                               ent(i,j,nz) * tr(i,j,nz-1))
       endif ; enddo
-      do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+      do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
         tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
       endif ; enddo ; enddo
     enddo
diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90
index 850480e3e6..3d2b322759 100644
--- a/src/tracer/MOM_tracer_hor_diff.F90
+++ b/src/tracer/MOM_tracer_hor_diff.F90
@@ -728,7 +728,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
   ! mixed and buffer layer corresponds to, such that:
   !     GV%Rlay(max_kRho-1) < Rml_max <= GV%Rlay(max_kRho)
   !$OMP parallel do default(shared) private(k_min,k_max,k_test)
-  do j=js-2,je+2 ; do i=is-2,ie+2 ; if (G%mask2dT(i,j) > 0.5) then
+  do j=js-2,je+2 ; do i=is-2,ie+2 ; if (G%mask2dT(i,j) > 0.0) then
     if ((Rml_max(i,j) > GV%Rlay(nz)) .or. (nkmb+1 > nz)) then ; max_kRho(i,j) = nz+1
     elseif ((Rml_max(i,j) <= GV%Rlay(nkmb+1)) .or. (nkmb+2 > nz)) then ; max_kRho(i,j) = nkmb+1
     else
@@ -756,7 +756,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
   !$OMP parallel default(shared) private(ns,tmp,itmp)
   !$OMP do
   do j=js-1,je+1
-    do k=1,nkmb ; do i=is-1,ie+1 ; if (G%mask2dT(i,j) > 0.5) then
+    do k=1,nkmb ; do i=is-1,ie+1 ; if (G%mask2dT(i,j) > 0.0) then
       if (h(i,j,k) > h_exclude) then
         num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j)
         k0_srt(i,ns,j) = k
@@ -764,7 +764,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
         h_srt(i,ns,j) = h(i,j,k)
       endif
     endif ; enddo ; enddo
-    do k=nkmb+1,PEmax_kRho ; do i=is-1,ie+1 ; if (G%mask2dT(i,j) > 0.5) then
+    do k=nkmb+1,PEmax_kRho ; do i=is-1,ie+1 ; if (G%mask2dT(i,j) > 0.0) then
       if ((k<=k_end_srt(i,j)) .and. (h(i,j,k) > h_exclude)) then
         num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j)
         k0_srt(i,ns,j) = k
@@ -812,7 +812,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
 !$OMP                                  kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, &
 !$OMP                                  wt_b,left_set,right_set,h_supply_frac_R,     &
 !$OMP                                  h_supply_frac_L)
-  do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.5) then
+  do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then
     ! Set up the pairings for fluxes through the zonal faces.
 
     do k=1,num_srt(i,j)   ; h_demand_L(k) = 0.0 ; h_used_L(k) = 0.0 ; enddo
@@ -965,7 +965,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
 !$OMP                                  kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, &
 !$OMP                                  wt_b,left_set,right_set,h_supply_frac_R,     &
 !$OMP                                  h_supply_frac_L)
-  do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.5) then
+  do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then
     ! Set up the pairings for fluxes through the meridional faces.
 
     do k=1,num_srt(i,j)   ; h_demand_L(k) = 0.0 ; h_used_L(k) = 0.0 ; enddo
@@ -1120,7 +1120,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
 !$OMP                          private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,Tr_La, &
 !$OMP                                     Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,h_L,h_R, &
 !$OMP                                     Tr_flux,Tr_adj_vert,wt_a,vol)
-      do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.5) then
+      do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then
         ! Determine the fluxes through the zonal faces.
 
         ! Find the acceptable range of tracer concentration around this face.
@@ -1260,7 +1260,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
 !$OMP                          private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,             &
 !$OMP                                  Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,        &
 !$OMP                                  h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol)
-      do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.5) then
+      do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then
         ! Determine the fluxes through the meridional faces.
 
         ! Find the acceptable range of tracer concentration around this face.
@@ -1377,7 +1377,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
 !$OMP                                  tr_flux_conv,Tr_flux_3d,k0a_Lv,Tr_adj_vert_L,&
 !$OMP                                  deep_wt_Rv,k0a_Rv,Tr_adj_vert_R) &
 !$OMP                          private(kLa,kLb,kRa,kRb,wt_b,wt_a)
-      do i=is,ie ; do J=js-1,je ; if (G%mask2dCv(i,J) > 0.5) then
+      do i=is,ie ; do J=js-1,je ; if (G%mask2dCv(i,J) > 0.0) then
         do k=1,nPv(i,J)
           kLb = k0b_Lv(J)%p(i,k); kRb = k0b_Rv(J)%p(i,k)
           if (deep_wt_Lv(J)%p(i,k) >= 1.0) then
@@ -1402,7 +1402,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
       endif ; enddo ; enddo
 !$OMP parallel do default(none) shared(PEmax_kRho,is,ie,js,je,G,h,Tr,tr_flux_conv,m)
       do k=1,PEmax_kRho ; do j=js,je ; do i=is,ie
-        if ((G%mask2dT(i,j) > 0.5) .and. (h(i,j,k) > 0.0)) then
+        if ((G%mask2dT(i,j) > 0.0) .and. (h(i,j,k) > 0.0)) then
           Tr(m)%t(i,j,k) = Tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / &
                                             (h(i,j,k)*G%areaT(i,j))
           tr_flux_conv(i,j,k) = 0.0
diff --git a/src/tracer/boundary_impulse_tracer.F90 b/src/tracer/boundary_impulse_tracer.F90
index f85623bbd1..4f5a3ea873 100644
--- a/src/tracer/boundary_impulse_tracer.F90
+++ b/src/tracer/boundary_impulse_tracer.F90
@@ -270,7 +270,7 @@ subroutine boundary_impulse_tracer_column_physics(h_old, h_new, ea, eb, fluxes,
 
   ! Set surface conditions
   do m=1,1
-    if (CS%remaining_source_time>0.0) then
+    if (CS%remaining_source_time > 0.0) then
       do k=1,CS%nkml ; do j=js,je ; do i=is,ie
         CS%tr(i,j,k,m) = 1.0
       enddo ; enddo ; enddo
diff --git a/src/tracer/ideal_age_example.F90 b/src/tracer/ideal_age_example.F90
index 5913251b14..a809e77a4c 100644
--- a/src/tracer/ideal_age_example.F90
+++ b/src/tracer/ideal_age_example.F90
@@ -352,7 +352,7 @@ subroutine ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G,
           exp((year-CS%tracer_start_year(m)) * CS%sfc_growth_rate(m))
     endif
     do k=1,CS%nkml ; do j=js,je ; do i=is,ie
-      if (G%mask2dT(i,j) > 0.5) then
+      if (G%mask2dT(i,j) > 0.0) then
         CS%tr(i,j,k,m) = sfc_val
       else
         CS%tr(i,j,k,m) = CS%land_val(m)