function p = psim(z)
% 	momentum aerodynamic profile stability correction
if z<0
	x1=ones(size(z),1);
	x=(1 - 15 *z).^0.25;
	p = log( ( x1+x.*x ) .* ( (x1+x).^2 ) / 8) - 2*atan(x) + pi/2;
else
	p=-4.7 * z;
end