function psi(x) REAL(rp), INTENT(IN) :: x REAL(rp) :: psi psi = 0.5_rp*(ERF(x) - 2.0_rp*x*EXP(-x**2)/SQRT(C_PI))/x**2 end function psi