;+
;-
;--------------------------------------------------------
function psit, zet
;--------------------------------------------------------

x=(1.-15*zet)^.5
psik=2*alog((1+x)/2)
x=(1.-34.15*zet)^.3333
psic=1.5*alog((1+x+x*x)/3)-sqrt(3.)*atan((1+2*x)/sqrt(3.))+4*atan(1.)/sqrt(3.)
f=zet*zet/(1.+zet*zet)
psi=(1.-f)*psik+f*psic  
   
ii=where(zet ge 0.)

if (ii(0) ne -1) then begin
  c=min([50,.35*zet])
  psi(ii)=-((1+2./3*zet(ii))^1.5+.6667*(zet(ii)-14.28)/exp(c(ii))+8.525)
endif

return, psi

end
;--------------------------------------------------------