c-----*----------------------------------------------------6---------7--
c     Correlation coefficient r between two series 
c       x(i) and y(i) with missing data, where i=1,...,n.
c     input: x(n),y(n),undef
c       undef is missing data; however, if ax or sx is still undefined,
c       it will be -9.99e33 whatever the input undef is.
c     output: r and nr
c       nr is the effective number of independent and defined values 
c       simultianeuosly in two series. 	
c     By Dr. LI Jianping, Decmber 5, 2001.
      subroutine correlationmiss(n,x,y,undef,r,nr)
      dimension x(n),y(n)
      dimension xw(n),yw(n)  ! work array
      undef1=-9.99e33
      num=0
      do 10 i=1,n
        if(x(i).eq.undef.or.y(i).eq.undef) goto 10
        num=num+1
        xw(num)=x(i)
        yw(num)=y(i)
  10  continue
      if(num.le.5)then
        r=undef1
      else
        call correlation(num,xw(1),yw(1),r)
      endif
      nr=num
      return
      end
c-----*----------------------------------------------------6---------7--
c     For the correlation coefficient r between two series 
c       x(i) and y(i), where i=1,...,n.
c     input: n,x(n),y(n)
c       n: number of time series
c       x(n): raw series 
c       y(n): raw series 
c     output: r
c       r: correlation coefficient between x and y
c     By Dr. LI Jianping, January 5, 2000.
      subroutine correlation(n,x,y,r)
      dimension x(n),y(n)
      call meanvar(n,x,ax,sx,vx)
      call meanvar(n,y,ay,sy,vy)
      sxy=0.
      do 10 i=1,n
        sxy=sxy+(x(i)-ax)*(y(i)-ay)
  10  continue
      sxy=sxy/float(n)
      r=sxy/(sx*sy)
      return
      end
c-----*----------------------------------------------------6---------7--
c     Computing the mean ax, standard deviation sx
c       and variance vx of a series x(i) (i=1,...,n).
c     input: n and x(n)
c       n: number of raw series
c       x(n): raw series
c     output: ax, sx and vx
c       ax: the mean value of x(n)
c       sx: the standard deviation of x(n)
c       vx: the variance of x(n)
c     By Dr. LI Jianping, May 6, 1998.
      subroutine meanvar(n,x,ax,sx,vx)
      dimension x(n)
      ax=0.
      vx=0.
      sx=0.
      do 10 i=1,n
        ax=ax+x(i)
  10  continue
      ax=ax/float(n)
      do 20 i=1,n
        vx=vx+(x(i)-ax)**2
  20  continue
      vx=vx/float(n)
      sx=sqrt(vx)
      return
      end