_RL htol ! stability parameter at zth [-]
htol = huol*ht/hu
xsq = SQRT( ABS(exf_one - htol*16. _d 0) )
xsq = MAX(SQRT(ABS(exf_one - 16.*htol)),exf_one)
psixh = -psim_fac*htol*stable + (exf_one-stable)*
xsq = MAX(SQRT(ABS(exf_one - htol*16. _d 0)),exf_one)
psixh = -psim_fac*htol*stable
htol = huol*zth/zwd
xsq = SQRT( ABS(one - htol*16. _d 0) )
psixh = -5. _d 0*htol*stable