!gslib gauinv.f
module gauinv_mod
implicit none
integer, parameter :: wp = selected_real_kind(10,50)
real(kind=wp),parameter :: lim=1.0e-10
real(kind=wp),parameter :: p0=-0.322232431088
real(kind=wp),parameter :: p1=-1.0
real(kind=wp),parameter :: p2=-0.342242088547
real(kind=wp),parameter :: p3=-0.0204231210245
real(kind=wp),parameter :: p4=-0.0000453642210148
real(kind=wp),parameter :: q0=0.0993484626060
real(kind=wp),parameter :: q1=0.588581570495
real(kind=wp),parameter :: q2=0.531103462366
real(kind=wp),parameter :: q3=0.103537752850
real(kind=wp),parameter :: q4=0.0038560700634

contains



      subroutine gauinv(p,xp,ierr)
!-----------------------------------------------------------------------
!
! Computes the inverse of the standard normal cumulative distribution
! function with a numerical approximation from : Statistical Computing,
! by W.J. Kennedy, Jr. and James E. Gentle, 1980, p. 95.
!
!
!
! INPUT/OUTPUT:
!
!   p    = double precision cumulative probability value: dble(psingle)
!   xp   = G^-1 (p) in single precision
!   ierr = 1 - then error situation (p out of range), 0 - OK
!
!
!-----------------------------------------------------------------------
     real(kind=wp) :: p
!
! Check for an error situation:
!
      ierr = 1
      if(p.lt.lim) then
            xp = -1.0e10
            return
      end if
      if(p.gt.(1.0-lim)) then
            xp =  1.0e10
            return
      end if
      ierr = 0      
!
! Get k for an error situation:
!
      pp   = p
      if(p.gt.0.5) pp = 1 - pp
      xp   = 0.0
      if(p.eq.0.5) return
!
! Approximate the function:
!
      y  = sqrt(log(1.0/(pp*pp)))
      xp = real( y + ((((y*p4+p3)*y+p2)*y+p1)*y+p0) / &
                    ((((y*q4+q3)*y+q2)*y+q1)*y+q0) )
      if(real(p).eq.real(pp)) xp = -xp
!
! Return with G^-1(p):
!
      return
      end subroutine
end module
