module kind_module
  implicit none
  save

  integer, parameter :: int_kind = KIND(1)
#ifdef unicos
  integer, parameter :: real_kind = KIND(1.0)
#else
  integer, parameter :: real_kind = KIND(1.0d0)
#endif

end module kind_module

module numbers_module
  use kind_module
  implicit none
  save

  private
  public :: zero, one, two, three, half

  real(real_kind), parameter :: zero=0.0, one=1.0, two=2.0, three=3.0, half=0.5

end module numbers_module

module clock_module
  use kind_module
  implicit none

  private
  public :: clock

contains

  subroutine clock (time)

    ! Out:
    real(real_kind), intent(out) :: time

    ! Local:
#ifdef cray
    real :: SECOND

    time = SECOND()
#else
    real :: user, syst, cputot, ETIME

    cputot = ETIME(user, syst)
    time = user + syst
#endif

    return
  end subroutine clock

end module clock_module

module type_module
  use kind_module
  implicit none

  private
  public :: pipe_type

  type pipe_type
     integer(int_kind) :: ncells
     real(real_kind), dimension(:), pointer :: pressure
     real(real_kind), dimension(:), pointer :: temp
     real(real_kind), dimension(:), pointer :: velocity
     real(real_kind), dimension(:), pointer :: density
     real(real_kind), dimension(:), pointer :: mass
  end type pipe_type

end module type_module

module data_module
  use kind_module

  implicit none
  save

  private
  public :: pressure, temp, density, velocity, mass

  real(real_kind), dimension(:), allocatable :: pressure, temp, density, velocity, mass
end module data_module

module comp_module
  use kind_module
  use type_module

  implicit none
  save

  private
  public :: comp, pipe_type

  type(pipe_type), dimension(:), allocatable :: comp
end module comp_module

module crunch1_module
  implicit none

  private
  public :: crunch1

contains

  subroutine crunch1 (pressure, temp, velocity, density, mass)
    use kind_module
    use numbers_module

    integer(int_kind) :: ncells
    real(real_kind), dimension(:), intent(in   ) :: velocity, density
    real(real_kind), dimension(:), intent(inout) :: pressure, temp
    real(real_kind), dimension(:), intent(  out) :: mass
    real(real_kind), dimension(:), allocatable :: vcell

    ncells = size(pressure)
    allocate (vcell(ncells))

    temp = temp*density
    pressure = pressure*temp
    vcell = velocity(1:ncells)
    vcell = half*(vcell + velocity(2:ncells+1))
    mass = density*vcell + pressure

    deallocate (vcell)

    return
  end subroutine crunch1
end module crunch1_module

module crunch2_module
  implicit none

contains

  subroutine crunch2 (jstart, jend)
    use kind_module
    use numbers_module
    use data_module

    integer(int_kind), intent(in) :: jstart, jend
    integer(int_kind) :: ncells
    real(real_kind), dimension(:), allocatable :: vcell

    ncells = jend - jstart + 1
    allocate (vcell(ncells))

    temp(jstart:jend) = temp(jstart:jend)*density(jstart:jend)
    pressure(jstart:jend) = pressure(jstart:jend)*temp(jstart:jend)
    vcell = velocity(jstart:jend)
    vcell = half*(vcell + velocity(jstart+1:jend+1))
    mass(jstart:jend) = density(jstart:jend)*vcell + pressure(jstart:jend)

    deallocate (vcell)

    return
  end subroutine crunch2
end module crunch2_module

module crunch3_module
  implicit none

  private
  public :: crunch3

contains

  subroutine crunch3 (pipe)
    use kind_module
    use type_module
    use numbers_module

    type(pipe_type), intent(inout) :: pipe
    integer(int_kind) :: ncells
    real(real_kind), dimension(:), allocatable :: vcell

    ncells = pipe%ncells
    allocate (vcell(ncells))

    pipe%temp = pipe%temp*pipe%density
    pipe%pressure = pipe%pressure*pipe%temp
    vcell = pipe%velocity(1:ncells)
    vcell = half*(vcell + pipe%velocity(2:ncells+1))
    pipe%mass = pipe%density*vcell + pipe%pressure

    deallocate (vcell)

    return
  end subroutine crunch3
end module crunch3_module

module crunch4_module
  implicit none

  private
  public :: crunch4

contains

  subroutine crunch4 (icomp)
    use kind_module
    use type_module
    use numbers_module
    use comp_module

    integer(int_kind), intent(in) :: icomp
    integer(int_kind) :: ncells
    real(real_kind), dimension(:), allocatable :: vcell

    ncells = comp(icomp)%ncells
    allocate (vcell(ncells))

    comp(icomp)%temp = comp(icomp)%temp*comp(icomp)%density
    comp(icomp)%pressure = comp(icomp)%pressure*comp(icomp)%temp
    vcell = comp(icomp)%velocity(1:ncells)
    vcell = half*(vcell + comp(icomp)%velocity(2:ncells+1))
    comp(icomp)%mass = comp(icomp)%density*vcell + comp(icomp)%pressure

    deallocate (vcell)

    return
  end subroutine crunch4
end module crunch4_module

program speed
  use kind_module
  use type_module
  use numbers_module
  use clock_module
  use crunch1_module
  use crunch2_module
  use crunch3_module
  use crunch4_module
  use data_module
  use comp_module

  implicit none

  integer(int_kind), parameter :: ncomp=7, repeat=1000
  real(real_kind), parameter :: factor=100000.0

  integer(int_kind) :: i, j, ntot
  real(real_kind) :: start, end, speedup
  integer(int_kind), dimension(:), allocatable :: jstart, jend
  real(real_kind), dimension(10) :: time
  real(real_kind), dimension(:), allocatable :: pressure_init, temp_init, density_init, velocity_init, mass_init
  real(real_kind), dimension(:), allocatable :: pressure_ref, temp_ref, density_ref, velocity_ref, mass_ref
  real(real_kind), dimension(:), allocatable :: pressure_diff, temp_diff, density_diff, velocity_diff, mass_diff

  allocate (comp(ncomp), jstart(ncomp), jend(ncomp))
  print *
  print *, 'number of components: ', ncomp
  ntot = 0
  jstart(1) = 1
  do i=1,ncomp
     comp(i)%ncells = sqrt(factor*i)
     allocate (comp(i)%pressure(comp(i)%ncells))
     allocate (comp(i)%temp(comp(i)%ncells))
     allocate (comp(i)%density(comp(i)%ncells))
     allocate (comp(i)%mass(comp(i)%ncells))
     allocate (comp(i)%velocity(comp(i)%ncells+1))
     ntot = ntot + comp(i)%ncells
     jend(i) = ntot
     if (i /= ncomp) jstart(i+1) = ntot+1
  end do

  do i=1,ncomp
     write(6,1) '  component ', i, ' ==> ncells: ', comp(i)%ncells, ' jstart: ', jstart(i), ' jend: ', jend(i)
  end do
  print *
  print *, 'total number of cells: ', ntot
  print *

  allocate (pressure(ntot), temp(ntot), density(ntot), mass(ntot), velocity(ntot+1))
  allocate (pressure_init(ntot), temp_init(ntot), density_init(ntot), mass_init(ntot), velocity_init(ntot+1))
  allocate (pressure_ref(ntot), temp_ref(ntot), density_ref(ntot), mass_ref(ntot), velocity_ref(ntot+1))
  allocate (pressure_diff(ntot), temp_diff(ntot), density_diff(ntot), mass_diff(ntot), velocity_diff(ntot+1))

  ! Initialize "properties" and use for all cases.
  call random_number (pressure_init)
  call random_number (temp_init)
  call random_number (density_init)
  call random_number (mass_init)
  call random_number (velocity_init)

  ! Global property arrays, sections passed via arg list.

  pressure = pressure_init
  temp = temp_init
  density = density_init
  mass = mass_init
  velocity = velocity_init

  call clock (start)
  do j=1,repeat
     do i=1,ncomp
        call crunch1 (pressure(jstart(i):jend(i)),  &
                      temp(jstart(i):jend(i)),  &
                      velocity(jstart(i):jend(i)+1),  &
                      density(jstart(i):jend(i)),  &
                      mass(jstart(i):jend(i)))
     end do
  end do
  call clock (end)
  time(1) = end - start
  write(6,2) 'global properties, passed via arg:     ', time(1)

  ! Consider this result to be a reference.  All runs should agree.
  pressure_ref = pressure
  temp_ref = temp
  density_ref = density
  mass_ref = mass
  velocity_ref = velocity

  ! Global property arrays in module.  Indices on which to work passed via arg list.

  pressure = pressure_init
  temp = temp_init
  density = density_init
  mass = mass_init
  velocity = velocity_init

  call clock (start)
  do j=1,repeat
     do i=1,ncomp
        call crunch2 (jstart(i), jend(i))
     end do
  end do
  call clock (end)
  time(2) = end - start
  speedup = time(1) / time(2)
  write(6,2) 'global properties, data in module:     ', time(2), ',  speedup: ', speedup

  pressure_diff = pressure - pressure_ref
  temp_diff = temp - temp_ref
  density_diff = density - density_ref
  mass_diff = mass - mass_ref
  velocity_diff = velocity - velocity_ref
  do i=1,ntot
     if (pressure_diff(i) /= zero) print *, 'pressure difference: ', i, pressure_diff(i)
     if (temp_diff(i) /= zero) print *, 'temp difference: ', i, temp_diff(i)
     if (density_diff(i) /= zero) print *, 'density difference: ', i, density_diff(i)
     if (mass_diff(i) /= zero) print *, 'mass difference: ', i, mass_diff(i)
     if (velocity_diff(i) /= zero) print *, 'velocity difference: ', i, velocity_diff(i)
  end do

  ! Properties pointered arrays as members of component derived types.  Entire
  ! component passed to work routine.

  do i=1,ncomp
     comp(i)%pressure = pressure_init(jstart(i):jend(i))
     comp(i)%temp = temp_init(jstart(i):jend(i))
     comp(i)%density = density_init(jstart(i):jend(i))
     comp(i)%mass = mass_init(jstart(i):jend(i))
     comp(i)%velocity = velocity_init(jstart(i):jend(i)+1)
  end do
  call clock (start)
  do j=1,repeat
     do i=1,ncomp
        call crunch3 (comp(i))
     end do
  end do
  call clock (end)
  time(3) = end - start
  speedup = time(1) / time(3)
  write(6,2) 'components passed via arg:             ', time(3), ',  speedup: ', speedup

  do i=1,ncomp
     pressure_diff(jstart(i):jend(i)) = comp(i)%pressure - pressure_ref(jstart(i):jend(i))
     temp_diff(jstart(i):jend(i)) = comp(i)%temp - temp_ref(jstart(i):jend(i))
     density_diff(jstart(i):jend(i)) = comp(i)%density - density_ref(jstart(i):jend(i))
     mass_diff(jstart(i):jend(i)) = comp(i)%mass - mass_ref(jstart(i):jend(i))
     velocity_diff(jstart(i):jend(i)+1) = comp(i)%velocity - velocity_ref(jstart(i):jend(i)+1)
  end do
  do i=1,ntot
     if (pressure_diff(i) /= zero) print *, 'pressure difference: ', i, pressure_diff(i)
     if (temp_diff(i) /= zero) print *, 'temp difference: ', i, temp_diff(i)
     if (density_diff(i) /= zero) print *, 'density difference: ', i, density_diff(i)
     if (mass_diff(i) /= zero) print *, 'mass difference: ', i, mass_diff(i)
     if (velocity_diff(i) /= zero) print *, 'velocity difference: ', i, velocity_diff(i)
  end do

  ! Properties pointered arrays as members of component derived types.  Properties
  ! passed via arg list.  Note that same work routine as first test is used.

  do i=1,ncomp
     comp(i)%pressure = pressure_init(jstart(i):jend(i))
     comp(i)%temp = temp_init(jstart(i):jend(i))
     comp(i)%density = density_init(jstart(i):jend(i))
     comp(i)%mass = mass_init(jstart(i):jend(i))
     comp(i)%velocity = velocity_init(jstart(i):jend(i)+1)
  end do
  call clock (start)
  do j=1,repeat
     do i=1,ncomp
        call crunch1 (comp(i)%pressure, comp(i)%temp, comp(i)%velocity, comp(i)%density, comp(i)%mass)
     end do
  end do
  call clock (end)
  time(4) = end - start
  speedup = time(1) / time(4)
  write(6,2) 'components, properties passed via arg: ', time(4), ',  speedup: ', speedup

  do i=1,ncomp
     pressure_diff(jstart(i):jend(i)) = comp(i)%pressure - pressure_ref(jstart(i):jend(i))
     temp_diff(jstart(i):jend(i)) = comp(i)%temp - temp_ref(jstart(i):jend(i))
     density_diff(jstart(i):jend(i)) = comp(i)%density - density_ref(jstart(i):jend(i))
     mass_diff(jstart(i):jend(i)) = comp(i)%mass - mass_ref(jstart(i):jend(i))
     velocity_diff(jstart(i):jend(i)+1) = comp(i)%velocity - velocity_ref(jstart(i):jend(i)+1)
  end do
  do i=1,ntot
     if (pressure_diff(i) /= zero) print *, 'pressure difference: ', i, pressure_diff(i)
     if (temp_diff(i) /= zero) print *, 'temp difference: ', i, temp_diff(i)
     if (density_diff(i) /= zero) print *, 'density difference: ', i, density_diff(i)
     if (mass_diff(i) /= zero) print *, 'mass difference: ', i, mass_diff(i)
     if (velocity_diff(i) /= zero) print *, 'velocity difference: ', i, velocity_diff(i)
  end do

  ! Properties pointered arrays as members of component derived types.
  ! Array of all components resides in module.

  do i=1,ncomp
     comp(i)%pressure = pressure_init(jstart(i):jend(i))
     comp(i)%temp = temp_init(jstart(i):jend(i))
     comp(i)%density = density_init(jstart(i):jend(i))
     comp(i)%mass = mass_init(jstart(i):jend(i))
     comp(i)%velocity = velocity_init(jstart(i):jend(i)+1)
  end do
  call clock (start)
  do j=1,repeat
     do i=1,ncomp
        call crunch4 (i)
     end do
  end do
  call clock (end)
  time(5) = end - start
  speedup = time(1) / time(5)
  write(6,2) 'components passed via arg:             ', time(5), ',  speedup: ', speedup

  do i=1,ncomp
     pressure_diff(jstart(i):jend(i)) = comp(i)%pressure - pressure_ref(jstart(i):jend(i))
     temp_diff(jstart(i):jend(i)) = comp(i)%temp - temp_ref(jstart(i):jend(i))
     density_diff(jstart(i):jend(i)) = comp(i)%density - density_ref(jstart(i):jend(i))
     mass_diff(jstart(i):jend(i)) = comp(i)%mass - mass_ref(jstart(i):jend(i))
     velocity_diff(jstart(i):jend(i)+1) = comp(i)%velocity - velocity_ref(jstart(i):jend(i)+1)
  end do
  do i=1,ntot
     if (pressure_diff(i) /= zero) print *, 'pressure difference: ', i, pressure_diff(i)
     if (temp_diff(i) /= zero) print *, 'temp difference: ', i, temp_diff(i)
     if (density_diff(i) /= zero) print *, 'density difference: ', i, density_diff(i)
     if (mass_diff(i) /= zero) print *, 'mass difference: ', i, mass_diff(i)
     if (velocity_diff(i) /= zero) print *, 'velocity difference: ', i, velocity_diff(i)
  end do

  stop
1 format (a, i2, a, i4, a, i5, a, i5)
2 format (a, f7.4, a, f6.3)
end program speed