!! Calculates the first three probability weighted moments of a sample,
!! using either the unbiased estimator (when A=B=0)
!! or a plotting position formula (when A<=B<=0).
!! This are alpha PWMs, following Rao & Hamed 2000, eqs. 3.1.4 and 3.1.6
Subroutine pwm (series, n, pwms, A, B, isBeta)
  integer :: n, isBeta
  real :: series (0:n-1), pwms(0:2), A, B
  integer :: i
  real :: acum(0:2), F

  acum (:) = 0. 
  if ( (A .eq. 0.) .and. (B .eq. 0.) ) then     !! use unbiased estimator
    do i = 1, n
      acum(0) = acum (0) + series (i-1)
      if (isBeta .eq. 0) then                   !! compute alpha PWMs
        acum(1) = acum(1) + series (i - 1) * (n - i)
        acum(2) = acum(2) + series (i - 1) * (n - i) * (n - i - 1)
      endif

      if (isBeta .eq. 1) then                   !! compute beta PWMs
        acum(1) = acum(1) + series (i - 1) * (i - 1)
        acum(2) = acum(2) + series (i - 1) * (i - 1) * (i - 2)
      endif
    enddo
  else  !! use plotting-position (biased) estimator
    do i = 1, n
      acum(0) = acum (0) + series (i-1)
      F = (real(i) + A) / (real(n) + B)
      if (isBeta .eq. 0) then                   !! compute alpha PWMs
        acum (1) = acum (1) + series(i - 1) * (1 - F)
        acum (2) = acum (2) + series(i - 1) * (1 - F) * (1 - F)
      endif
      if (isBeta == 1) then                     !! compute beta PWMs
        acum (1) = acum (1) + series (i - 1) * F
        acum (2) = acum (2) + series (i - 1) * F * F
      endif
    enddo
  endif

  pwms (0) = acum (0) / n
  pwms (1) = acum (1) / n / (n - 1)
  pwms (2) = acum (2) / n / ((n-1)*(n-2))

  return    
End Subroutine pwm

!Subroutine lMoments (series, n, lMoment, A, B)
!  integer :: n
!  real :: series (n), lMoment (n), A, B
!  integer :: i, j , orderLMom
!  real :: C, D, E, acum(3), alpha (3)
!
!  Call pwm (series,n, alpha, A, B, 0)
!
!  lMoment (2) = alpha (1)
!  lMoment (3) = alpha (1) - 2*alpha (2)
!  lMoment (4) = alpha (1) - 6*alpha (2)
!
!  return
!End Subroutine lMoments