c ----------------------------------------------------------------------
      subroutine escros (x,y,npts,xmsg,ymsg,xmean,ymean,xsd,ysd
     1                  ,mxlag,ccv,ccr,ier)
 
 
c this routine will estimate the crosscovariances and crosscorrelations
c .   for the input series x and y. missing data is allowed.
 
c .   to get the (+) and (-) ccv and ccr two calls must be made. the
c .   first with the calling sequence above and the second with the
c .   x and y series reversed (also the other arguments) . note
c .   that  the second ccv(0) and ccr(0) are redundant.
 
c arguments :
c .   x,y      - input vectors with which the calculations are to take place.
c .   npts     - length of vectors x and y
c .   xmsg,ymsg- missing code: if no msg values set to some number
c .                            which will not be encountered.
c .              xmsg will be used to fill missing values
c .   xmean    - mean of the x vector (if xmean=xmsg: program will calculate)
c .   ymean    - mean of the y vector (if ymean=ymsg: program will calculate)
c .   xsd      - st. dev. of x vector (if xsd  =xmsg: program will calculate)
c .   ysd      - st. dev. of y vector (if ysd  =ymsg: program will calculate)
c .   mxlag    - max lag crosscovariance/correlation to be estimated
c .   ccv      - vector of length .ge. (mxlag+1) containing sample cross cov
c .   ccr      - vector of length .ge. (mxlag+1) containing sample cross cor
 
c .                ccv(0)     = xyvar (lag 0  )  : ccr(0)     =  ( lag 0   )
c .                ccv(1)     =       (lag 1  )  : ccr(1)     =  ( lag 1   )
c .                           .                               .
c .                           .                               .
c .                ccv(mxlag) =       (lag mxlag): ccr(mxlag) =  ( lag mxlag)
 
c .   ier      - if (ier.ne.0) an error has occurred
 
c note : uncalculated quantities are set to xmsg
 
      double precision x(npts) , y(npts) , ccv(0:mxlag) , ccr(0:mxlag)
      double precision xmean, xvar, xsd, xmsg, ymean, yvar, ysd, ymsg
      double precision xyn, xy1
 
      ier = 0
      if (npts.lt.2) ier = 1
      if (mxlag.lt.0 .or. mxlag.gt.npts) ier = 2
      if ((xsd.ne.xmsg .and. xsd.le.0.d0) .or.
     *    (ysd.ne.ymsg .and. ysd.le.0.d0)) ier = 3
      do lag=0,max0( 0, min0(mxlag,npts) )
         ccv(lag) = xmsg
         ccr(lag) = xmsg
      enddo
      if (ier.ne.0) return
 
      if (xmean.eq.xmsg .or. xsd.eq.xmsg) then
         call stat2 (x,npts,xmsg,xmean,xvar,xsd,nptused,jer)
         if (jer.ne.0) then
             ier = -jer
         elseif (xsd.eq.0.d0) then
             ier = -5              ! x must be a series of constant values
         elseif (nptused.eq.0) then
             ier = -6              ! x must be a series of missing values
         endif
         if (ier.ne.0) return
      endif
 
      if (ymean.eq.ymsg .or. ysd.eq.ymsg) then
         call stat2 (y,npts,ymsg,ymean,yvar,ysd,nptused,jer)
         if (jer.ne.0) then
             ier = -(jer+100)
         elseif (ysd.eq.0.d0) then
             ier = -105            ! y must be a series of constant values
         elseif (nptused.eq.0) then
             ier = -106            ! y must be a series of missing values
         endif
         if (ier.ne.0) return
      endif
 
      xsdysd = 1./(xsd*ysd)
 
      do lag=0,mxlag
         xyn = 0.d0
         xy1 = 0.d0
       do n=1,npts-lag
          if (x(n).ne.xmsg .and. y(n+lag).ne.ymsg) then
              xy1 = xy1 + (y(n+lag)-ymean)*(x(n)-xmean)
              xyn = xyn + 1.d0
          endif
       enddo
         if (xyn.ge.2.) then
             ccv(lag) = xy1/(xyn-1.d0)
             ccr(lag) = ccv(lag)*xsdysd
         else
             ccv(lag) = xmsg
             ccr(lag) = xmsg
         endif
      enddo
 
      return
      end
c ----------------------------------------------------------------------
      subroutine stat2 (x,npts,xmsg,xmean,xvar,xsd,nptused,ier)
 
c this routine will calculate estimates of the first two  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
c .   xsd      - sqrt( xvar )
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)
 
      integer npts  , nptused, ier
      double precision x(npts), xmsg, xmean, xvar, xsd
      double precision x1, x2, x3, xn
 
      xmean   = xmsg
      xvar    = xmsg
      xsd     = xmsg
      nptused = 0
 
      if (npts.lt.1) then
         ier = 1
         return
      endif
 
      ier = 0
      x1  = 0.d0
      x2  = 0.d0
      xn  = 0.d0
      do n=1,npts
         if (x(n).ne.xmsg) then
             xn = xn + 1.0d0
             x1 = x1 + x(n)
             x2 = x2 + x(n)*x(n)
         endif
      enddo
 
      nptused = xn
 
      if ((xn-1.d0).gt.0.d0) then
         x3    = x1*x1/xn
c c c    xvar  = amax1( (x2-x3)/(xn-1.d0) , 0.d0 ) ! prevent possible roundoff
         xsd   = sqrt( xvar )
         xmean = x1/xn
      elseif ((xn-1.d0).eq.0.d0) then
         xmean = x1
         xvar  = 0.d0
         xsd   = 0.d0
      else
         ier = 2    ! error code for all msg values
      endif
 
      return
      end

