subroutine VariabilityIndicies(nsamp,xvar,nthresh,threshold,
     &shannon,coefvar,xmean,total,unevenness,standunevenness,
     &arearatio,pvindex,meaninterval,maxinterval,xintensity)
      implicit none
      integer nsamp,i,nthresh,j,k
      double precision xvar(nsamp),shannon,coefvar,xmean,total,
     &unevenness,standunevenness,arearatio,cumline(nsamp),
     &evenline(nsamp),areacum,areaeven,trans(nsamp),avetrans,
     &pvindex,meaninterval(nthresh),maxinterval(nthresh),
     &finterval(nsamp),fstart,threshold(nthresh),recumline(nsamp),
     &repvi,retrans(nsamp),averetrans,xintensity
!
      cumline(1)=xvar(1)
      recumline(1)=xvar(nsamp)
      k=0
      if(xvar(1).gt.0.0d0)k=1
      do i=2,nsamp
        cumline(i)=cumline(i-1)+xvar(i)
        recumline(i)=recumline(i-1)+xvar(nsamp-i+1)
        if(xvar(i).gt.0.0d0)k=k+1
      enddo
      total=cumline(nsamp)
      xmean=total/dble(nsamp)
      if(k.gt.0)then
        xintensity=total/dble(k)
      else
        xintensity=-9999.0d0
      endif
      if(nsamp.eq.1.or.total.le.0.0d0)then
        shannon=-9999.0d0
        coefvar=-9999.0d0
        unevenness=-9999.0d0
        standunevenness=-9999.0d0
        arearatio=-9999.0d0
        pvindex=-9999.0d0
        repvi=-9999.0d0
        do i=1,nthresh
          meaninterval(i)=-9999.0d0
          maxinterval(i)=-9999.0d0
        enddo
        return
      endif
      do i=1,nsamp
        evenline(i)=xmean*dble(i)
      enddo
      unevenness=0.0d0
      shannon=0.0d0
      coefvar=0.0d0
      avetrans=0.0d0
      averetrans=0.0d0
      do i=1,nsamp
        unevenness=unevenness+(cumline(i)-evenline(i))*
     &(cumline(i)-evenline(i))
        trans(i)=(cumline(i)-evenline(i))/evenline(i)
!        trans(i)=2.0d0*(cumline(i)-evenline(i))/(xmean*dble((nsamp+1)))
        retrans(i)=(recumline(i)-evenline(i))/evenline(i)
        avetrans=avetrans+trans(i)
        averetrans=averetrans+retrans(i)
        if(dabs(xvar(i)).gt.1.0d-10)then
          shannon=shannon+(xvar(i)/total)*dlog(xvar(i)/total)
        endif
        coefvar=coefvar+(xvar(i)-xmean)*(xvar(i)-xmean)
      enddo

!      avetrans=0.0d0
!      do i=1,nsamp
!        trans(i)=trans(i)*retrans(nsamp-i+1)
!        if(trans(i).lt.0.0d0)then
!          trans(i)=-dsqrt(-trans(i))
!        else
!          trans(i)=dsqrt(trans(i))
!        endif
!        avetrans=avetrans+trans(i)
!      enddo

      coefvar=dsqrt(coefvar/dble(nsamp-1))/xmean
      unevenness=dsqrt(unevenness/dble(nsamp))
      standunevenness=unevenness/xmean
      shannon=-shannon/dlog(dble(nsamp))
      avetrans=avetrans/dble(nsamp)
      averetrans=averetrans/dble(nsamp)
      pvindex=0.0d0
      repvi=0.0d0
      do i=1,nsamp
        pvindex=pvindex+(trans(i)-avetrans)*(trans(i)-avetrans)
        repvi=repvi+(retrans(i)-averetrans)*(retrans(i)-averetrans)
      enddo
      pvindex=dsqrt(pvindex/dble(nsamp))
      repvi=dsqrt(repvi/dble(nsamp))

!      pvindex=(pvindex+repvi)/2.0d0

      areacum=0.0d0
      areaeven=0.0d0
      do i=2,nsamp
        areacum=areacum+(cumline(i)+cumline(i-1))*0.5d0
        areaeven=areaeven+(evenline(i)+evenline(i-1))*0.5d0
      enddo
      arearatio=areacum/areaeven
!intervals
      do i=1,nthresh
        fstart=0.0d0
        k=1
        do j=1,nsamp
          if(xvar(j).ge.threshold(i))then
            finterval(k)=dble(j)-fstart-1.0d0
            fstart=dble(j)
            k=k+1
          endif
        enddo
        if((dble(nsamp)-fstart).gt.0.5d0)then
          finterval(k)=dble(nsamp)-fstart
        else
          k=k-1
        endif
        maxinterval(i)=finterval(1)
        meaninterval(i)=finterval(1)
        do j=2,k
          if(finterval(j).gt.maxinterval(i))maxinterval(i)=finterval(j)
          meaninterval(i)=meaninterval(i)+finterval(j)
        enddo
        meaninterval(i)=meaninterval(i)/dble(k)
      enddo
      return
      end subroutine VariabilityIndicies