diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90
index 6b60769144..051ac446db 100644
--- a/src/tracer/MOM_neutral_diffusion.F90
+++ b/src/tracer/MOM_neutral_diffusion.F90
@@ -1574,34 +1574,39 @@ real function interpolate_for_nondim_position(dRhoNeg, Pneg, dRhoPos, Ppos)
 
   character(len=120) :: mesg
 
-  if (Ppos < Pneg) then
-    call MOM_error(FATAL, 'interpolate_for_nondim_position: Houston, we have a problem! Ppos<Pneg')
-  elseif (dRhoNeg>dRhoPos) then
-    write(stderr,*) 'dRhoNeg, Pneg, dRhoPos, Ppos=',dRhoNeg, Pneg, dRhoPos, Ppos
-    write(mesg,*) 'dRhoNeg, Pneg, dRhoPos, Ppos=', dRhoNeg, Pneg, dRhoPos, Ppos
-    call MOM_error(WARNING, 'interpolate_for_nondim_position: '//trim(mesg))
-  elseif (dRhoNeg>dRhoPos) then !### Does this duplicated test belong here?
-    call MOM_error(FATAL, 'interpolate_for_nondim_position: Houston, we have a problem! dRhoNeg>dRhoPos')
-  endif
-  if (Ppos<=Pneg) then ! Handle vanished or inverted layers
-    interpolate_for_nondim_position = 0.5
-  elseif ( dRhoPos - dRhoNeg > 0. ) then
-    interpolate_for_nondim_position = min( 1., max( 0., -dRhoNeg / ( dRhoPos - dRhoNeg ) ) )
-  elseif ( dRhoPos - dRhoNeg == 0) then
-    if (dRhoNeg>0.) then
-      interpolate_for_nondim_position = 0.
-    elseif (dRhoNeg<0.) then
-      interpolate_for_nondim_position = 1.
-    else ! dRhoPos = dRhoNeg = 0
+  if ((Ppos > Pneg) .and. (dRhoPos - dRhoNeg >= 0. )) then
+    if ( dRhoPos - dRhoNeg > 0. ) then
+      interpolate_for_nondim_position = min( 1., max( 0., -dRhoNeg / ( dRhoPos - dRhoNeg ) ) )
+    elseif (dRhoPos - dRhoNeg == 0) then
+      if (dRhoNeg > 0.) then
+        interpolate_for_nondim_position = 0.
+      elseif (dRhoNeg < 0.) then
+        interpolate_for_nondim_position = 1.
+      else ! dRhoPos = dRhoNeg = 0
+        interpolate_for_nondim_position = 0.5
+      endif
+    else ! dRhoPos - dRhoNeg < 0
       interpolate_for_nondim_position = 0.5
     endif
-  else ! dRhoPos - dRhoNeg < 0
+  elseif (Ppos == Pneg) then ! Handle vanished or inverted layers
     interpolate_for_nondim_position = 0.5
+  else ! ((Ppos < Pneg) .or. (dRhoNeg > dRhoPos) )
+    ! Error handling for problematic cases.  It is expected that this should never occur.
+    write(mesg,*) 'dRhoNeg, Pneg, dRhoPos, Ppos', dRhoNeg, Pneg, dRhoPos, Ppos
+    call MOM_error(WARNING, 'interpolate_for_nondim_position: '//trim(mesg))
+    !  write(stderr,*) trim(mesg)
+    if ((Ppos < Pneg) .and. (dRhoNeg > dRhoPos)) then
+      mesg = '(Ppos < Pneg) and (dRhoNeg > dRhoPos)'
+    elseif (Ppos < Pneg) then
+      mesg = 'Ppos < Pneg'
+    elseif (dRhoNeg > dRhoPos) then
+      mesg = trim(mesg)//'; dRhoNeg > dRhoPos'
+    else  ! This should never happen.
+      mesg = 'Unexpected failure.'
+    endif
+    call MOM_error(FATAL, 'interpolate_for_nondim_position: '//trim(mesg))
   endif
-  if ( interpolate_for_nondim_position < 0. ) &
-    call MOM_error(FATAL, 'interpolate_for_nondim_position: Houston, we have a problem! Pint < Pneg')
-  if ( interpolate_for_nondim_position > 1. ) &
-    call MOM_error(FATAL, 'interpolate_for_nondim_position: Houston, we have a problem! Pint > Ppos')
+
 end function interpolate_for_nondim_position
 
 !> Higher order version of find_neutral_surface_positions. Returns positions within left/right columns