function Psi = Psi_m(zoL)
unstable = find(zoL<0);
stable   = find(zoL>0);
neutral  = find(zoL==0);

Psi(neutral) = 0;

x = (1-16*zoL).^.25;
Psi(unstable) = 2*log( (1+x(unstable))/2 ) + log( (1+x(unstable).^2)/2 ) - 2*atan(x(unstable)) + pi/2 ; 

Psi(stable) = -5*zoL(stable);