SUBROUTINE nscor1 (s, n, n2, work, ifault)
 
! Code converted using TO_F90 by Alan Miller
! Date: 2003-06-17  Time: 12:13:20

!     Algorithm AS 177   Appl. Statist. (1982) Vol. 31, No. 2

!     Exact calculation of Normal Scores

IMPLICIT NONE
INTEGER, PARAMETER   :: dp = SELECTED_REAL_KIND(15, 100)

INTEGER, INTENT(IN)        :: n2
REAL (dp), INTENT(OUT)     :: s(n2)
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: work(4,721)
INTEGER, INTENT(OUT)       :: ifault

INTERFACE
  SUBROUTINE init(work)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(15, 100)
    REAL (dp), INTENT(OUT)  :: work(4,721)
  END SUBROUTINE init
END INTERFACE

REAL (dp) :: c, scor, ai1, ani, an

INTEGER   :: i, i1, j, ni
REAL (dp), PARAMETER  :: one = 1.0D0, zero = 0.0D0, h = 0.025D0
INTEGER, PARAMETER    :: nstep = 721

ifault=3
IF (n2 /= n/2) RETURN
ifault=1
IF (n <= 1) RETURN
ifault=0
IF (n > 2000) ifault=2

an=n
c=LOG(an)

!        Accumulate ordinates for calculation of integral for rankits

DO  i=1, n2
  i1=i-1
  ni=n-i
  ai1=i1
  ani=ni
  scor=zero
  DO  j=1,nstep
    scor = scor + EXP(work(2,j) + ai1*work(3,j) + ani*work(4,j) + c) * work(1,j)
  END DO
  s(i)=scor * h
  c=c + LOG(ani/i)
END DO
RETURN

END SUBROUTINE nscor1



SUBROUTINE init(work)

!        Algorithm AS 177.1   Appl. Statist. (1982) Vol. 31, No. 2

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(15, 100)

INTERFACE
  FUNCTION alnorm( x, upper ) RESULT( fn_val )
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(15, 100)
    REAL(dp), INTENT(IN)   ::  x
    LOGICAL, INTENT(IN)    ::  upper
    REAL(dp)               ::  fn_val
  END FUNCTION alnorm
END INTERFACE

REAL (dp), INTENT(OUT)  :: work(4,721)

REAL (dp)  :: xx
REAL (dp), PARAMETER  :: xstart = -9.0D0, h = 0.025D0, pi2 = -0.918938533D0, &
                         half = 0.5D0
INTEGER, PARAMETER    :: nstep = 721
INTEGER  :: i

xx=xstart

!        Set up arrays for calculation of integral

DO  i=1,nstep
  work(1,i)=xx
  work(2,i)=pi2 - xx * xx * half
  work(3,i)=LOG(alnorm(xx, .true.))
  work(4,i)=LOG(alnorm(xx, .false.))
  xx=xstart + i * h
END DO
RETURN
END SUBROUTINE init



!  Algorithm AS66 Applied Statistics (1973) vol.22, no.3

!  Evaluates the tail area of the standardised normal curve
!  from x to infinity if upper is .true. or
!  from minus infinity to x if upper is .false.

! ELF90-compatible version by Alan Miller
! Latest revision - 29 November 2001

FUNCTION alnorm( x, upper ) RESULT( fn_val )
IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(15, 100)

REAL(DP), INTENT(IN)   ::  x
LOGICAL, INTENT(IN)    ::  upper
REAL(DP)               ::  fn_val

!  Local variables
REAL(DP), PARAMETER   ::  zero=0.0_DP, one=1.0_DP, half=0.5_DP, con=1.28_DP
REAL(DP)              ::  z, y
LOGICAL               ::  up

!  Machine dependent constants
REAL(DP), PARAMETER  ::  ltone = 7.0_DP, utzero = 18.66_DP
REAL(DP), PARAMETER  ::  p = 0.398942280444_DP, q = 0.39990348504_DP,   &
                         r = 0.398942280385_DP, a1 = 5.75885480458_DP,  &
                         a2 = 2.62433121679_DP, a3 = 5.92885724438_DP,  &
                         b1 = -29.8213557807_DP, b2 = 48.6959930692_DP, &
                         c1 = -3.8052E-8_DP, c2 = 3.98064794E-4_DP,     &
                         c3 = -0.151679116635_DP, c4 = 4.8385912808_DP, &
                         c5 = 0.742380924027_DP, c6 = 3.99019417011_DP, &
                         d1 = 1.00000615302_DP, d2 = 1.98615381364_DP,  &
                         d3 = 5.29330324926_DP, d4 = -15.1508972451_DP, &
                         d5 = 30.789933034_DP

up = upper
z = x
IF( z < zero ) THEN
   up = .NOT. up
   z = -z
END IF
IF( z <= ltone .OR. (up .AND. z <= utzero) ) THEN
   y = half*z*z
   IF( z > con ) THEN
      fn_val = r*EXP( -y )/(z+c1+d1/(z+c2+d2/(z+c3+d3/(z+c4+d4/(z+c5+d5/(z+c6))))))
   ELSE
      fn_val = half - z*(p-q*y/(y+a1+b1/(y+a2+b2/(y+a3))))
   END IF
ELSE
   fn_val = zero
END IF

IF( .NOT. up ) fn_val = one - fn_val
RETURN
END FUNCTION alnorm



!----------------------------------------------------------------------------


SUBROUTINE nscor2(s, n, n2, ier)

!     Algorithm as 177.3, applied statistics, v.31, 161-165, 1982.

!     Calculates approximate expected values of normal order statistics.
!     claimed accuracy is 0.0001, though usually accurate to 5-6 dec.

!     Arguments:
!     s(n2)   = output, the first n2 expected values.
!     n       = input, the sample size.
!     n2      = input, the number of order statistics required; must
!                      be <= n/2.
!     ier     = output, error indicator
!                   = 0 if no error detected
!                   = 1 if n <= 1.
!                   = 2 if n > 2000, in which case the order statistics
!                          are still calculated, but may be inaccurate.
!                   = 3 if n2 > n/2 (n.b. this differs from the
!                          published algorithm which returns an error
!                          if n2 is not equal to n/2.)

!     Calls ppnd = Applied Statistics algorithm 111.
!     An alternative is ppnd7 in algorithm AS 241.

! N.B. This routine is not in double precision.

IMPLICIT NONE

INTEGER, INTENT(IN)   :: n2
REAL, INTENT(OUT)     :: s(n2)
INTEGER, INTENT(IN)   :: n
INTEGER, INTENT(OUT)  :: ier

INTERFACE
  FUNCTION correc(i, n) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, INTENT(IN)  :: i
    INTEGER, INTENT(IN)  :: n
    REAL                 :: fn_val
  END FUNCTION correc

  SUBROUTINE ppnd7 (p, normal_dev, ifault)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: sp = SELECTED_REAL_KIND(6, 30)
    REAL (sp), INTENT(IN)   :: p
    INTEGER, INTENT(OUT)    :: ifault
    REAL (sp), INTENT(OUT)  :: normal_dev
  END SUBROUTINE ppnd7
END INTERFACE

REAL     :: an, ai, e1, e2, l1
INTEGER  :: i, k
REAL , PARAMETER  ::   &
    eps(4) = (/ 0.419885, 0.450536, 0.456936, 0.468488 /),  &
    dl1(4) = (/ 0.112063, 0.121770, 0.239299, 0.215159 /),  &
    dl2(4) = (/ 0.080122, 0.111348, -0.211867, -0.115049 /),  &
    gam(4) = (/ 0.474798, 0.469051, 0.208597, 0.259784 /),  &
    lam(4) = (/ 0.282765, 0.304856, 0.407708, 0.414093 /),  &
    bb = -0.283833, d = -0.106136, b1 = 0.5641896

!     input parameter checks.

ier = 3
IF(n2 > n/2) RETURN
ier = 1
IF(n <= 1) RETURN
ier = 0
IF(n > 2000) ier = 2
s(1) = b1
IF(n == 2) RETURN

!     calculate normal tail areas for first 3 order statistics.

an = n
k = 3
IF(n2 < k) k = n2
DO  i = 1,k
  ai = i
  e1 = (ai - eps(i))/(an + gam(i))
  e2 = e1**lam(i)
  s(i) = e1 + e2*(dl1(i) + e2*dl2(i))/an - correc(i,n)
END DO
IF(n2 == k) GO TO 20

!     calculate normal areas for other cases.

DO  i = 4,n2
  ai = i
  l1 = lam(4) + bb/(ai + d)
  e1 = (ai - eps(4))/(an + gam(4))
  e2 = e1**l1
  s(i) = e1 + e2*(dl1(4) + e2*dl2(4))/an - correc(i,n)
END DO

!     convert tail areas to normal deviates.

20 DO  i = 1,n2
  CALL ppnd7(s(i), s(i), ier)
END DO

RETURN
END SUBROUTINE nscor2



FUNCTION correc(i, n) RESULT(fn_val)

! Calculates correction for tail area of the i-th largest of n order statistics.

IMPLICIT NONE

INTEGER, INTENT(IN)  :: i
INTEGER, INTENT(IN)  :: n
REAL                 :: fn_val

REAL  :: an
REAL, PARAMETER  :: c1(7) = (/ 9.5, 28.7, 1.9, 0., -7.0, -6.2, -1.6 /),  &
    c2(7) = (/ -6195., -9569., -6728., -17614., -8278., -3570., 1075. /),  &
    c3(7) = (/ 9.338E4, 1.7516E5, 4.1040E5, 2.1576E6, 2.376E6, 2.065E6,  &
    2.065E6 /), mic = 1.e-6, c14 = 1.9E-5

fn_val = c14
IF(i*n == 4) RETURN
fn_val = 0.0
IF(i < 1 .OR. i > 7) RETURN
IF(i /= 4 .AND. n > 20) RETURN
IF(i == 4 .AND. n > 40) RETURN
an = n
an = 1.0/(an*an)
fn_val = (c1(i) + an*(c2(i) + an*c3(i)))*mic
RETURN
END FUNCTION correc



SUBROUTINE ppnd7 (p, normal_dev, ifault)

! ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3, 477- 484.

! Produces the normal deviate Z corresponding to a given lower tail area of P;
! Z is accurate to about 1 part in 10**7.

! The hash sums below are the sums of the mantissas of the coefficients.
! They are included for use in checking transcription.

! This ELF90-compatible version by Alan Miller - 20 August 1996
! N.B. The original algorithm is as a function; this is a subroutine

IMPLICIT NONE
INTEGER, PARAMETER  :: sp = SELECTED_REAL_KIND(6, 30)

REAL (sp), INTENT(IN)   :: p
INTEGER, INTENT(OUT)    :: ifault
REAL (sp), INTENT(OUT)  :: normal_dev

! Local variables

REAL (sp) :: zero = 0.0, one = 1.0, half = 0.5, split1 = 0.425,  &
             split2 = 5.0, const1 = 0.180625, const2 = 1.6, q, r

! Coefficients for P close to 0.5

REAL (sp) :: a0 = 3.3871327179E+00, a1 = 5.0434271938E+01, &
                   a2 = 1.5929113202E+02, a3 = 5.9109374720E+01, &
                   b1 = 1.7895169469E+01, b2 = 7.8757757664E+01, &
                   b3 = 6.7187563600E+01
! HASH SUM AB          32.3184577772

! Coefficients for P not close to 0, 0.5 or 1.

REAL (sp) :: c0 = 1.4234372777E+00, c1 = 2.7568153900E+00, &
             c2 = 1.3067284816E+00, c3 = 1.7023821103E-01, &
             d1 = 7.3700164250E-01, d2 = 1.2021132975E-01
! HASH SUM CD    15.7614929821

! Coefficients for P near 0 or 1.

REAL (sp) :: e0 = 6.6579051150E+00, e1 = 3.0812263860E+00, &
             e2 = 4.2868294337E-01, e3 = 1.7337203997E-02, &
             f1 = 2.4197894225E-01, f2 = 1.2258202635E-02
! HASH SUM EF    19.4052910204

ifault = 0
q = p - half
IF (ABS(q) <= split1) THEN
  r = const1 - q * q
  normal_dev = q * (((a3 * r + a2) * r + a1) * r + a0) / &
               (((b3 * r + b2) * r + b1) * r + one)
  RETURN
ELSE
  IF (q < zero) THEN
    r = p
  ELSE
    r = one - p
  END IF
  IF (r <= zero) THEN
    ifault = 1
    normal_dev = zero
    RETURN
  END IF
  r = SQRT(-LOG(r))
  IF (r <= split2) THEN
    r = r - const2
    normal_dev = (((c3 * r + c2) * r + c1) * r + c0) / ((d2 * r + d1) * r + one)
  ELSE
    r = r - split2
    normal_dev = (((e3 * r + e2) * r + e1) * r + e0) / ((f2 * r + f1) * r + one)
  END IF
  IF (q < zero) normal_dev = - normal_dev
  RETURN
END IF
RETURN
END SUBROUTINE ppnd7