! example for specifying the precision for program, functions, subroutines
! at only one place. Declare variables as real(pr), specify all real constants
! such as 0.1 in the form 0.1_pr

module prec
  implicit none
!     Change the following line to change precision everywhere
  integer, parameter :: pr=kind(1.) 
!     choices for specifying pr:
!     (1) values 4, 8, 16 on DEC Alpha for single, double, quad precision.
!         disadvantage: program might not work on other machine/compiler
!     (2) pr = kind(1.)  for default single precision   
!         pr = kind(1.D0) for default double precision
!     (3) pr = selected_real_kind(7,99)
!         gives lowest precision numbers which have at least 7 digits of
!         decimal precision, and at least a range from 1.E-99 to 1.E99
end module

module vectorfun
  use prec
  implicit none
contains
  function length(v)        ! computes the euclidian length of a vector
    implicit none
    real(pr)    :: length, v(:), s
    integer :: n, i

    n = size(v)
    s = 0._pr
    do i=1,n
      s = s + v(i)**2
    end do
    length=sqrt(s)
  end function
end module

program test
  use vectorfun 
  use prec
  implicit none

  integer :: i,n
  real(pr)    :: r 
  real(pr), allocatable :: x(:)

! print *,'base =',radix(r),' digits =',digits(r),' emin =',minexponent(r), &
!      ' emax =',maxexponent(r),' epsilon (chopping) =',epsilon(r), &
!      ' largest =',huge(r),' smallest =',tiny(r) 

  print *,'n ='
  read *, n
  allocate(x(n))
  do i=1,n
    print *,'x(',i,') ='
    read *, x(i)
  end do
  r=length(x)
  print *, 'length =',r
end program