!************************************************************
!*  This program calculates the statistical moments of a    *
!*  distribution: Mean, Variance, Skewness, etc.            *
!* -------------------------------------------------------- *
!* Ref.: "Numerical Recipes, by W.H. Press, B.P. Flannery,  *
!*        S.A. Teukolsky and T. Vetterling, Cambridge       *
!*        University Press, 1986".                          *
!*                                                          *
!*                    F90 Release By J-P Moreau, Paris.     *
!* -------------------------------------------------------- *
!* SAMPLE RUN:                                              *
!*                                                          *
!* Number of data: 5                                        *
!*   1: 12                                                  *
!*   2: 9                                                   *
!*   3: 7                                                   *
!*   4: 15                                                  *
!*   5: 6                                                   *
!*                                                          *
!* Average ..........:   9.80000000000000                   *
!* Average  Deviation:   2.96000000000000                   *
!* Standard Deviation:   3.70135110466435                   *
!* Variance .........:   13.7000000000000                   *
!* Skewness .........:  0.291548695888740                   *
!* Kurtosis .........:  -1.90779903031595                   *
!*                                                          *  
!************************************************************
Program TMOMENT

  integer ndata
  real*8, pointer ::  Y(:)
  real*8 average,avdev,stdev,variance,skewness,kurtosis
  
  print *,' '
  write(*,10,advance='no'); read *, ndata
  allocate(Y(1:ndata),stat=ialloc)
  do i=1, ndata
    write(*,20,advance='no') i
	read *, Y(i)
  end do

  call Moment(Y,ndata,average,avdev,stdev,variance,skewness,kurtosis)

! print results
  print *,' '
  print *,' Average ..........: ', average
  print *,' Average  Deviation: ', avdev
  print *,' Standard Deviation: ', stdev
  print *,' Variance .........: ', variance
  print *,' Skewness .........: ', skewness
  print *,' Kurtosis .........: ', kurtosis
  print *,' '
  stop

10 format('  Number of data: ')
20 format(' ',i2,': ')

End


!************************************************************
!* Given an array of Data of length n, this routine returns * 
!* its mean ave, average deviation adev, standard deviation *
!* sdev, variance var, skewness skew, and kurtosis curt.    *
!************************************************************
Subroutine Moment(data,n,ave,adev,sdev,var,skew,curt)
  integer n
  real*8 data(n), ave,adev,sdev,var,skew,curt, p,s
  if(n.le.1) then
    print *,' N must be at least 2!'
	stop
  end if
  s=0.d0
  do j=1,n
    s=s+data(j)
  end do
! calculate mean
  ave=s/n
  adev=0.d0
  var=0.d0
  skew=0.d0
  curt=0.d0
  do j=1,n
    s=data(j)-ave
	adev=adev+dabs(s)
	p=s*s
	var=var+p
	p=p*s
	skew=skew+p
	p=p*s
	curt=curt+p
  end do
  adev=adev/n
  var=var/(n-1)
  sdev=dsqrt(var)
  if(var.ne.0.d0) then
    skew=skew/(n*sdev**3)
	curt=curt/(n*var**2)-3.d0
  else
    print *,' No skew or kurtosis when zero variance.'
  end if
  return
end

!end of file tmoment.f90