MODULE lmdz_arth

  IMPLICIT NONE; PRIVATE
  PUBLIC arth

  INTEGER, PARAMETER :: NPAR_ARTH = 16, NPAR2_ARTH = 8

  INTERFACE arth
    ! Returns an arithmetic progression, given a first term "first", an
    ! increment and a number of terms "n" (including "first").

    MODULE PROCEDURE arth_r, arth_i
    ! The difference between the procedures is the kind and type of
    ! arguments first and increment and of function result.
  END INTERFACE

CONTAINS

  PURE FUNCTION arth_r(first, increment, n)

    REAL, INTENT(IN) :: first, increment
    INTEGER, INTENT(IN) :: n
    REAL arth_r(n)

    ! Local:
    INTEGER :: k, k2
    REAL :: temp

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

    IF (n > 0) arth_r(1) = first
    IF (n <= NPAR_ARTH) THEN
      DO k = 2, n
        arth_r(k) = arth_r(k - 1) + increment
      END DO
    ELSE
      DO k = 2, NPAR2_ARTH
        arth_r(k) = arth_r(k - 1) + increment
      END DO
      temp = increment * NPAR2_ARTH
      k = NPAR2_ARTH
      DO
        IF (k >= n) EXIT
        k2 = k + k
        arth_r(k + 1:min(k2, n)) = temp + arth_r(1:min(k, n - k))
        temp = temp + temp
        k = k2
      END DO
    END IF

  END FUNCTION arth_r

  !*************************************

  PURE FUNCTION arth_i(first, increment, n)

    INTEGER, INTENT(IN) :: first, increment, n
    INTEGER arth_i(n)

    ! Local:
    INTEGER :: k, k2, temp

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

    IF (n > 0) arth_i(1) = first
    IF (n <= NPAR_ARTH) THEN
      DO k = 2, n
        arth_i(k) = arth_i(k - 1) + increment
      END DO
    ELSE
      DO k = 2, NPAR2_ARTH
        arth_i(k) = arth_i(k - 1) + increment
      END DO
      temp = increment * NPAR2_ARTH
      k = NPAR2_ARTH
      DO
        IF (k >= n) EXIT
        k2 = k + k
        arth_i(k + 1:min(k2, n)) = temp + arth_i(1:min(k, n - k))
        temp = temp + temp
        k = k2
      END DO
    END IF

  END FUNCTION arth_i

END MODULE lmdz_arth