Update Thermohaline Mode Properties Solver

Sometimes we encounter failures to converge leading to this error message:

Lines 195 to 199 in 1baa434

    
           if((myvars(2)<0).or.(ierr /= 0)) then 
        
              write(*,*) "WARNING: thermohaline Newton relaxation failed to converge, falling back to estimate" 
        
              maxl2 = maxl*maxl 
        
           else ! NR succeeded, so use results in myvars 
        
              !Plug solution into "l^2" and lambda.

Adrian may have a solution to some of the underlying issues.

I claim that my python code here appropriately handles the if myvars(2)<0 case (what's happening is it it's trying to solve an equation that is quadratic in k^2 and therefore equivalently a quartic in k; if it solve the quadratic and finds k^2 < 0, then as a fallback it goes back and solves the full quartic to ensure k is real).

Rich suggests we might want to instead just do a bracketed search where we solve the full eigenvalue problem (only a 3x3 matrix fortunately) at every k and directly search for the k that maximizes the growth rate, rather than doing this (formally equivalent, but somewhat different in terms of numerics) root-finding business.

Rich has provided what we need! See 27a18c6 and c7d832a. We just need to point MESA's Brown and Harrington routines to these new solvers and test them.

After some testing, @rhdtownsend's solution seems to be working wonderfully. It is much more elegant and robust than the old 2D root find, so I think we can confidently mark this one as complete. Rich is writing up some notes to make sure we have this thoroughly documented.

	if((myvars(2)<0).or.(ierr /= 0)) then
	write(,) "WARNING: thermohaline Newton relaxation failed to converge, falling back to estimate"
	maxl2 = maxl*maxl
	else ! NR succeeded, so use results in myvars
	!Plug solution into "l^2" and lambda.