c ----------------------------------------------------------------------
      subroutine stat4 (x,npts,xmsg,xmean,xvar,xsd,xskew,xkurt
     1                 ,nptused,ier)
 
c this routine will calculate estimates of the first four moments
c .   of the vector x containing missing data.
 
c input arguments:
c .   x        - input vector
c .   npts     - length of x
c .   xmsg     - missing code: if there are no msg values
c .                            set xmsg to some value which will
c .                            not be encountered.
 
c output arguments:
c .   xmean    - mean of x (first moment)
c .   xvar     - sample variance of x (second moment)
c .   xsd      - sqrt( xvar )
c .   xskew    - coef. of skewness (third moment)
c .              departure from symmetry. if skew>0 [skew<0] the
c .              distribution trails off to the right [left] .
c .   xkurt    - coef. of kurtosis (fourth moment)
c .              the normal distribution has a kurtosis of 3 . this
c .              value is subtracted from the calculated kurtosis.
c .              thus neg values are possible and the returned value
c .              is kurtosis relative to the normal distribution.
c .              if (kurt-3) > 0 [<0] it is usually more sharply peaked
c .              [flatter] than the normal distribution. (leptokurtic and
c .              platykurtic , respectively.) e.g. , a rectangular function has
c .              a kurtosis of 1.8 or 1.8-3 = -1.2 relative to the normal
c .              distribution.
c .   nptused  - no. of points used to calculate the estimates
c .   ier      - if (ier.ne.0) an error has occurred
 
c note :
c .   uncalculated quantities are set to xmsg (nptused set to zero)
 
      real x(1:npts)
      double precision x1, x2, x3, x4, dxvar
 
      xmean = xmsg
      xvar  = xmsg
      xsd   = xmsg
      xskew = xmsg
      xkurt = xmsg
      nptused = 0
 
 
      if (npts.lt.1) then
         ier = 1
         return
      endif
 
      ier = 0
      x1  = 0.d0
      x2  = 0.d0
      xn  = 0.
      do n=1,npts
         if (x(n).ne.xmsg) then
             xn = xn + 1.0
             x1 = x1 + dble(x(n))
             x2 = x2 + dprod( x(n),x(n) )
         endif
      enddo
 
      nptused = xn
 
      if ((xn-1.).gt.0.) then
         x3    = x1*x1/dble(xn)
         dxvar = dmax1( (x2-x3)/dble(xn-1.) , 0.d0 ) !prevent possible roundoff
         xvar  = sngl( dxvar )
         xsd   = sqrt( xvar )
         xmean = sngl(x1/dble(xn))
 
         if (dxvar.gt.0.d0) then
             x1 = 0.d0
             x2 = 0.d0
             x3 = 0.d0
             x4 = 0.d0
             do n=1,npts
                if (x(n).ne.xmsg) then
                    x4 = dble(x(n)-xmean)
                    x1 = x4**3
                    x2 = x2 + x1
                    x3 = x3 + x1*x4
                endif
             enddo
             xskew = sngl( (x2/dsqrt(dxvar)**3)/dble(xn))
             xkurt = sngl( (x3/(dxvar*dxvar))/dble(xn)) - 3.
         endif
      elseif ((xn-1.).eq.0.) then
         xmean = sngl(x1)
         xvar  = 0.0
         xsd   = 0.0
      else
         ier = 2    ! error code for all msg values
         return
      endif
 
      return
      end

