WRITE(*, '(a, i4, a, i4)')' Replicate: ', replicate,  &
                              '  Percentile: ', percentile

    sumsq_LS = zero

    REWIND(unit_data)
    IF (line1 > 1) THEN
      DO i = 1, line1-1
        READ(unit_data, *)
      END DO
    END IF
    case = 1
    weight = minus1
    x(0) = one
    DO
      IF (ypos > nvar) THEN
        READ(unit_data, *, IOSTAT=iostatus) x(1:nvar), y
      ELSE IF (ypos == 1) THEN
        READ(unit_data, *, IOSTAT=iostatus) y, x(1:nvar)
      ELSE
        READ(unit_data, *, IOSTAT=iostatus) x(1:ypos-1), y, x(ypos:nvar)
      END IF

      IF (iostatus > 0) CYCLE                   ! Error in data
      IF (iostatus < 0) EXIT                    ! End of file

      IF(ANY(case == order(i1:i2))) THEN        ! Delete case if in this 10%
        CALL includ(weight, x, y)
      END IF
      case = case + 1
    END DO

! N.B. Subroutine INCLUD increases nobs even when weights are negative.
! Calculate correct value for the present nobs.

    nobs = nobs_full - (i2 + 1 - i1)

! Find subsets which fit well

    IF (fit_const) THEN
      nvar_max = max_size - 1
    ELSE
      nvar_max = max_size
    END IF
    CALL init_subsets(nvar_max, fit_const)
                             ! Re-order the QR-factorization if variables are
                             ! to be forced in or out.
    IF (first > 1) THEN
      CALL reordr(init_vorder, first-1, 1, ier)
    END IF
    IF (last < ncol) THEN
      CALL reordr(init_vorder, last, 1, ier)
    END IF
    SELECT CASE(search_method)
      CASE(1)
        CALL efroym(first, last, fin, fout, nsize, ier, lout)
      CASE(2)
        CALL seqrep(first, last, ier)
      CASE(3)
        CALL seq2(first, last, ier)
      CASE(4)
        CALL seq2(first, last, ier)
        CALL xhaust(first, last, ier)
    END SELECT

! Pick best subset

    min_crit_val = HUGE(one)
    IF (search_method > 1) THEN
      IF (criterion == 3) variance = sserr / (nobs - ncol)
      nsize = 0
      DO i = first, max_size
        IF (criterion /= 3) variance = ress(i,1) / (nobs - i)
        CALL calc_penalty(i, first, variance, criterion, penalty)
        crit_val = ress(i,1) + penalty
        IF (nsize == 0 .OR. crit_val < min_crit_val) THEN
          min_crit_val = crit_val
          nsize = i
        END IF
      END DO
    END IF

    ipos = ((nsize-1)*nsize)/2 + 1
    CALL reordr(lopt(ipos:,1), nsize, 1, ier)

! Estimate the regression coefficients using the LS-projections

    ALLOCATE( beta_LS(1:nsize), vorder_cpy(1:nsize) )
    vorder_cpy = vorder(1:nsize)
    CALL shell(vorder_cpy, nsize)                 ! Shell sort from find_sub
    WRITE(unit_rpt, 970) percentile, vorder_cpy(1:nsize)
    970 FORMAT('Percentile no.', i3, '   Selected variables:'/(' ', 15i5))
    CALL regcf(beta_LS, nsize, ier)
!    WRITE(unit_rpt, 980) beta_LS
!    980 FORMAT('LS regression coefficients:'/(' ', 6g13.5))
    vorder_cpy = vorder(1:nsize)

! Return the order of variables to the original order, as in the data set

    DO i = 1, nvar
      CALL reordr(list, i, 2, ier)
    END DO

! Estimate the 10% omitted, and re-instate the deleted cases

    REWIND(unit_data)
    IF (line1 > 1) THEN
      DO i = 1, line1-1
        READ(unit_data, *)
      END DO
    END IF
    case = 1
    weight = one
    x(0) = one
    DO
      IF (ypos > nvar) THEN
        READ(unit_data, *, IOSTAT=iostatus) x(1:nvar), y
      ELSE IF (ypos == 1) THEN
        READ(unit_data, *, IOSTAT=iostatus) y, x(1:nvar)
      ELSE
        READ(unit_data, *, IOSTAT=iostatus) x(1:ypos-1), y, x(ypos:nvar)
      END IF

      IF (iostatus > 0) CYCLE                   ! Error in data
      IF (iostatus < 0) EXIT                    ! End of file

      IF(ANY(case == order(i1:i2))) THEN        ! Restore case if in this 10%
        fit_LS = zero
        DO i = 1, nsize
          j = vorder_cpy(i)
          fit_LS = fit_LS + beta_LS(i) * x(j)
        END DO
        sumsq_LS = sumsq_LS + (y - fit_LS)**2
        CALL includ(weight, x, y)                ! INCLUD destroys x
      END IF
      case = case + 1
    END DO
    WRITE(unit_rpt,1000) sumsq_LS
    1000 FORMAT('Sums of sq. (LS) = ', g12.4)
    total_LS = total_LS + sumsq_LS

!   CALL print_QR
    DEALLOCATE( beta_LS, vorder_cpy )

  END DO             ! percentile = 1, 10

  WRITE(unit_rpt, 1030) total_LS
  WRITE(*, 1030) total_LS
  1030 FORMAT(/' Total sum of squares (LS) = ', g13.5)
  WRITE(unit_rpt, '(/)')
  msep = msep + total_LS

END DO             ! replicate = 1, nrepl

msep = msep / (nrepl * nobs_full)
WRITE(*, 900) msep
WRITE(unit_rpt, 900) msep
900 FORMAT(' Overall mean squared error of prediction = ', g13.5)
WRITE(*, 910) SQRT(msep)
WRITE(unit_rpt, 910) SQRT(msep)
910 FORMAT(' RMS (prediction error) = ', g13.5)

DEALLOCATE( order, x, list, seeds )
STOP
END SUBROUTINE cross_validation



SUBROUTINE ranord(order, n)

! Generate a random ordering of the integers 1 ... n.

INTEGER, INTENT(IN)  :: n
INTEGER, INTENT(OUT) :: order(:)

! Local variables

REAL      :: wk(n)
INTEGER   :: i

DO i = 1, n
  order(i) = i
END DO
CALL RANDOM_NUMBER(wk)

CALL sqsort(wk, n, order)

RETURN
END SUBROUTINE ranord




SUBROUTINE sqsort(a, n, t)

!   NON-RECURSIVE STACK VERSION OF QUICKSORT FROM N.WIRTH'S PASCAL
!   BOOK, 'ALGORITHMS + DATA STRUCTURES = PROGRAMS'.

!   SINGLE PRECISION, ALSO CHANGES THE ORDER OF THE ASSOCIATED ARRAY T.

INTEGER, INTENT(IN)     :: n
INTEGER, INTENT(IN OUT) :: t(:)
REAL, INTENT(IN OUT)    :: a(:)

! Local variables

INTEGER :: i, j, k, l, r, s, stackl(15), stackr(15), ww
REAL    :: w, x

s = 1
stackl(1) = 1

! KEEP TAKING THE TOP REQUEST FROM THE STACK UNTIL S = 0.

stackr(1) = n

10 l = stackl(s)
r = stackr(s)

! KEEP SPLITTING A(L), ... , A(R) UNTIL L>=R.

s = s - 1

20 i = l
j = r
k = (l + r)/2

! REPEAT UNTIL I > J.

x = a(k)

30 IF(a(i) >= x) GO TO 40
i = i + 1
GO TO 30

40 IF(x >= a(j)) GO TO 50
j = j - 1
GO TO 40

50 IF(i > j) GO TO 60
w = a(i)
ww = t(i)
a(i) = a(j)
t(i) = t(j)
a(j) = w
t(j) = ww
i = i + 1
j = j - 1
IF(i <= j) GO TO 30

60 IF(j - l < r - i) GO TO 75
IF(l >= j) GO TO 65
s = s + 1
stackl(s) = l
stackr(s) = j

65 l = i
GO TO 90

75 IF(i >= r) GO TO 80
s = s + 1
stackl(s) = i
stackr(s) = r

80 r = j

90 IF(l < r) GO TO 20

IF(s /= 0) GO TO 10

RETURN
END SUBROUTINE sqsort



SUBROUTINE calc_penalty(size1, size2, var, penalty_num, penalty_val)

! Calculate value of penalty for size of subset.
! Currently the penalties available are:
! Number Name
!    1   AIC
!    2   BIC
!    3   Mallows' Cp
!    4   Hannan-Quinn
!    5   Efroymson (F-to-delete = F-to-add = 4.0)

INTEGER, INTENT(IN)    :: size1, size2, penalty_num
REAL(dp), INTENT(IN)   :: var
REAL(dp), INTENT(OUT)  :: penalty_val

! Local variables
REAL(dp)  :: zero = 0.0_dp, one = 1.0_dp, two = 2.0_dp, four = 4.0_dp

IF (size1 == size2) THEN
  penalty_val = zero
  RETURN
END IF

IF (penalty_num < 1 .OR. penalty_num > 5) THEN
  penalty_val = zero
  RETURN
END IF

SELECT CASE(penalty_num)
  CASE(1)
    penalty_val = rss(size1) * (one - exp(two * (size2 - size1) / nobs))
  CASE(2)
    penalty_val = rss(size1) *    &
                  (one - exp((size2 - size1) * LOG(REAL(nobs)) / nobs))
  CASE(3)
    penalty_val = two * var * (size1 - size2)
  CASE(4)
    penalty_val = rss(size1) *    &
                  (one - exp(two * (size2 - size1) * LOG(LOG(REAL(nobs))) / nobs))
  CASE(5)
    penalty_val = four * var * (size1 - size2)
END SELECT

RETURN

END SUBROUTINE calc_penalty



SUBROUTINE get_numbers(list, n)
!     Read in a list of numbers which may be separated by commas, blanks
!     or either `..' or `-'.

INTEGER, DIMENSION(:), INTENT(OUT) :: list
INTEGER, INTENT(OUT)               :: n

!     Local variables
CHARACTER (LEN=4)   :: delimiters = ' ,-.'
CHARACTER (LEN=100) :: text
INTEGER             :: nmax, length, i1, i2, iostatus, i, number
LOGICAL             :: sequence

nmax = SIZE(list)

start: DO
  WRITE(*, *) 'Enter variable numbers on one line'
  WRITE(*, *) 'e.g. 1-5 8 11 .. 15  use commas or blanks as separators'
  WRITE(*, *) ': '
  READ(*, '(a)') text
  text = ADJUSTL(text)
  length = LEN_TRIM(text)
  IF (length == 0) THEN
    n = 0
    RETURN
  END IF

  n = 1
  i1 = 1
  sequence = .FALSE.
  DO
    i2 = SCAN( text(i1:), delimiters )
    IF (i2 == 0) THEN
      i2 = length
    ELSE
      i2 = i2 + i1 - 2
    END IF
    READ(text(i1:i2), *, IOSTAT=iostatus) number
    IF (iostatus /= 0) THEN
      WRITE(*, *) '** Error: numeric data expected **'
      WRITE(*, '(1x, a)') text(1:length)
      text = ' '
      DO i = i1, i2
        text(i:i) = '^'
      END DO
      WRITE(*, '(1x, a)') text(1:i2)
      CYCLE start
    END IF

    IF (sequence) THEN
      IF (number <= list(n-1)) THEN
        WRITE(*, *) 'Variable numbers not increasing'
        WRITE(*, '(1x, a)') text(1:length)
        text = ' '
        DO i = i1, i2
          text(i:i) = '^'
        END DO
        WRITE(*, '(1x, a)') text(1:i2)
        CYCLE start
      END IF
      DO
        list(n) = list(n-1) + 1
        IF (list(n) >= number) EXIT
        n = n + 1
      END DO
    ELSE
      list(n) = number
    END IF
    IF (i2 == length) RETURN

    i1 = i2 + 1
    sequence = .FALSE.

                                       ! Find end of delimiters
    DO
      IF ( SCAN( text(i1:i1), delimiters ) > 0) THEN
        IF (text(i1:i1) == '-' .OR. text(i1:i1+1) == '..') sequence = .TRUE.
        i1 = i1 + 1
      ELSE
        EXIT
      END IF
    END DO
    n = n + 1
    IF (n > nmax) THEN
      WRITE(*, *) '** Too many numbers entered - list truncated **'
      n = nmax
      RETURN
    END IF
  END DO
END DO start

RETURN
END SUBROUTINE get_numbers



SUBROUTINE set_seed()

INTEGER, ALLOCATABLE :: seed(:)
INTEGER              :: k

!     Set the random number seed.

CALL RANDOM_SEED(size=k)
ALLOCATE (seed(k))
CALL RANDOM_SEED(get=seed)
WRITE(*, *)'Old random number seeds: ', seed

WRITE(*, '(1x, a, i4, a)') 'Enter ', k, ' integers as random number seeds: '
READ(*, *) seed
WRITE(11, '(a/ 10(" ", i10))') 'New random number seeds:', seed
CALL RANDOM_SEED(put=seed)

RETURN
END SUBROUTINE set_seed

END PROGRAM subset