Better error messages in interp_for_nondim_pos

  Cleaned up the error messages when there are fatal problems in
interpolate_for_nondim_position() in the MOM_neutral_diffusion module to
explicitly give an indications of all the problems.  These had previously
triggered fatal errors (or should have), but had less explicit or incorrect
error messages in some cases, so all solutions and behavior are identical in any
cases that worked previously.

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