54 real ::
psih,x,z,zeta,l
55 real,
parameter :: a=1.,b=0.667,c=5.,d=0.35,eps=1.e-20
57 if ((l.ge.0).and.(l.lt.eps))
then
59 else if ((l.lt.0).and.(l.gt.(-1.*eps)))
then
63 if ((log10(z)-log10(abs(l))).lt.log10(eps))
then
68 psih = - (1.+0.667*a*zeta)**(1.5) - b*(zeta-c/d)*exp(-d*zeta) &
71 x=(1.-16.*zeta)**(.25)
72 psih=2.*log((1.+x*x)/2.)