c-----*----------------------------------------------------6---------7--
c     nt-month (or nt-season, or others) lagged and leading correlation coefficients
c        rt(-nt:nt,nm) between two anomaly series x(ny,nm) and y(ny,nm).
c     input: ny,nt,nm,x(ny,nm),y(ny,nm)
c       ny: year
c       nm: month, or season or others
c       nt: maximum lag or lead time (in month or season, etc.)
c       x(ny,nm) and y(ny,nm): raw anomaly series
c     output: rt(-nt:nt,nm)
c       rt(i,j): nt-month (or nt-season, or others) lagged and leading 
c              correlation coefficients between x and y
c              i<0, lagged correlations, x trails y;
c              i=0, simultaneous correltion bewteen x and y;
c              i>0, leaing correlations, x lead y.
c     By Dr. LI Jianping, June 8, 2005.
      subroutine mllcorrelation(ny,nm,nt,x,y,rt)
      dimension x(ny,nm),y(ny,nm)
      dimension rt(-nt:nt,nm)
c     Simultaneous correltion bewteen x and y.
      do j=1,nm
        call correlation(ny,x(1,j),y(1,j),rt(0,j))
      enddo
c     Leading correlations, x leads y.
      do j=1,nm
        do j1=1,nt
          j2=int((j1+j-1)/float(nm))+1
          j3=(j1+j)-nm*(j2-1)
          nj=ny-j2+1
c         write(100,*)j,j1,j2,j3,nj
          call correlation(nj,x(1,j),y(j2,j3),rt(j1,j))
        enddo
      enddo
c     Lagged correlations, x trails y.
      do j=nm,1,-1
        do j1=1,nt
          j2=int((j1+nm-j)/float(nm))+1
          j3=nm*j2-(j1+nm-j)
          nj=ny-j2+1
c         write(100,*)j,j1,j2,j3,nj
          call correlation(nj,x(j2,j),y(1,j3),rt(-j1,j))
        enddo
      enddo
      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