c----*----------------------------------------------------6---------7--
c     Correlations between two anomaly series of
c        12 calendar months, 4 seasons (DJF, MAM, JJA and SON),
c        monthly, seasonal and annual data from monthly anomaly series x(i,12)
c        and y(i,12) (i=1,...,n)
c     input: n, x(n,12), y(n,12)
c       n: year
c       x(n,12),y(n,12): the raw anomaly series
c     output: r(24)
c       r(j): correlation cofficients between x and y, j=1,...,24
c             j=1,12 for 12 calenar months, 
c             j=13,16 for 4 seasons (DJF, MAM, JJA, and SON)
c             j=17,20 for D-F, M-M, J-J, S-N month to month data
c             j=21,24 for montly, annual, NDJFMA,  and MJJASO data
c     By Dr. LI Jianping, February 6, 2004.
      subroutine mscorrelation_2(n,x,y,r)
      dimension x(n,12),y(n,12),r(24)
      call mcorrelation_1(n,x,y,r(1))   !1-12
      call sscorrelation_1(n,x,y,r(13)) !13-16
      call mmcorrelation_1(n,x,y,12,2,r(17))
      call mmcorrelation_1(n,x,y,3,5,r(18))
      call mmcorrelation_1(n,x,y,6,8,r(19))
      call mmcorrelation_1(n,x,y,9,11,r(20))
      call mmcorrelation_1(n,x,y,1,12,r(21))
      call scorrelation_1(n,x,y,1,12,r(22))
      call mmcorrelation_1(n,x,y,11,4,r(23))
      call mmcorrelation_1(n,x,y,5,10,r(24))
      return
      end
c-----*----------------------------------------------------6---------7--
c     For seasonal variation of the correlation coefficients r
c       between two monthly anomaly series x(i,12) and y(i,12), i=1,...,n.
c     input: n, x(n,12), y(n,12)
c       n: year
c       x(n,12),y(n,12): the raw anomaly series
c     output: r(13)
c       r(j): correlation cofficients between x and y, j=1,...,13
c             j=1,12 for 12 calenar months, 
c     By Dr. LI Jianping, September 6, 2001.
      subroutine mcorrelation_1(n,x,y,r)
      dimension x(n,12),y(n,12),r(13)
      do i=1,12
        call correlation(n,x(1,i),y(1,i),r(i))
      enddo
      r(13)=r(1)
      return
      end
c-----*----------------------------------------------------6---------7--
c     For the correlation coefficients r between two anomaly series
c       x(i,12) and y(i,12) for 4 seasons (DJF, MAM, JJA and SON), i=1,...,n.
c       The raw series are monthly.
c     input: n, x(n,12), y(n,12)
c       n: year
c       x(n,12),y(n,12): the raw anomaly series
c     output: r(4)
c       r(j): correlation cofficients between x and y for 4 seasons
c     By Dr. LI Jianping, Decmber 5, 2001.
      subroutine sscorrelation_1(n,x,y,r)
      dimension x(n),y(n),r(4)
      dimension zx1(n),zx2(n+1) !work array
      dimension zy1(n),zy2(n+1) !work array
      dimension mon0(4),mon1(4) !work array
      data mon0/12,3,6,9/
      data mon1/ 2,5,8,11/
      do i=1,4
        call sseries_1(n,x,mon0(i),mon1(i),zx1,zx2)
        call sseries_1(n,y,mon0(i),mon1(i),zy1,zy2)
        if(mon0(i).le.mon1(i))then
          call correlation(n,zx1,zy1,r(i))
        else
          call correlation(n+1,zx2,zy2,r(i))
        endif
      enddo
      return
      end
c-----*----------------------------------------------------6---------7--
c     Correlation between two month-month series of a specified season
c       The specified season is from mon0 through mon1.
c     input: n, x(n,12), y(n,12), mon0, mon1
c       n: year 
c       x(n,12),y(n,12): the raw anomaly series
c       mon0: initial month
c       mon1: end month
c     output: r
c       r: correlation cofficient between x and y 
c     By Dr. LI Jianping, September 6, 2001.
      subroutine mmcorrelation_1(n,x,y,mon0,mon1,r)
      dimension x(n,12),y(n,12)
      dimension xw(n*12),yw(n*12)  ! work array
      call mmseries_1(n,x,mon0,mon1,nn,xw)
      call mmseries_1(n,y,mon0,mon1,nn,yw)
      call correlation(nn,xw,yw,r)
      return
      end
c-----*----------------------------------------------------6---------7--
c     For the correlation coefficient r between two anomaly series
c       x(i,12) and y(i,12) of a specified season or year, i=1,...,n.
c       The raw series are monthly.
c     input: n, x(n,12), y(n,12), mon0 and mon1
c       n: year
c       x(n,12),y(n,12): the raw anomaly series
c       mon0: initial month
c       mon1: end month
c     output: r
      subroutine scorrelation_1(n,x,y,mon0,mon1,r)
      dimension x(n,12),y(n,12)
      dimension zx1(n),zx2(n+1) !work array
      dimension zy1(n),zy2(n+1) !work array
      call sseries_1(n,x,mon0,mon1,zx1,zx2)
      call sseries_1(n,y,mon0,mon1,zy1,zy2)
      if(mon0.le.mon1)then
        call correlation(n,zx1,zy1,r)
      else
        call correlation(n+1,zx2,zy2,r)
      endif
      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, 1999.
      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
c-----*----------------------------------------------------6---------7--
c     Reconstructing the series of a specified season or year from 
c       a monthly series x(i,12) (i=1,...,n).
c       The specified season is from mon0 through mon1.
c     input: n, x(n,12), mon0 and mon1
c       n: year
c       x(n,12): raw monthly series
c       mon0: initial month
c       mon1: end month
c     output: zx1 and zx2
c     The specified seaonal mean here is the average from mon0 to mon1. 
c     By Dr. LI Jianping, September 6, 2001.
      subroutine sseries_1(n,x,mon0,mon1,zx1,zx2)
      dimension x(n,12),zx1(n),zx2(n+1)
      do k=1,n
        zx1(k)=0.
        zx2(k)=0.
      enddo
      zx2(n+1)=0.
      if(mon0.le.mon1)then
        nk=mon1-mon0+1
        do k=1,n
          do i=mon0,mon1
            zx1(k)=zx1(k)+x(k,i) 
          enddo
          zx1(k)=zx1(k)/nk
          zx2(k)=zx1(k)
        enddo
      else
        nk=12+mon1-mon0+1
        do i=1,mon1
          zx2(1)=zx2(1)+x(1,i)
        enddo
        zx2(1)=zx2(1)/mon1
        do k=2,n
          do i=mon0,12
            zx2(k)=zx2(k)+x(k-1,i)
          enddo
          do i=1,mon1
            zx2(k)=zx2(k)+x(k,i)
          enddo
          zx2(k)=zx2(k)/nk
        enddo
        do i=mon0,12
          zx2(n+1)=zx2(n+1)+x(n,i)
        enddo
        zx2(n+1)=zx2(n+1)/(12-mon0+1)
        do k=1,n
          zx1(k)=zx2(k)
        enddo
      endif
      return
      end
c-----*----------------------------------------------------6---------7--
c     Computing the month-month series of a specified season from 
c       a series x(i,12) (i=1,...,n).
c       The specified season is from mon0 through mon1.
c     input: n, x(n,12), mon0 and mon1
c       n: year
c       x(n,12): raw monthly series
c       mon0: initial month
c       mon1: end month
c     output: nn,sx(nn)
c       nn: number of the series sx
c       sx: the monthly series for the specified season
c     By Dr. LI Jianping, September 6, 2001.
      subroutine mmseries_1(n,x,mon0,mon1,nn,sx)
      dimension x(n,12),sx(n*12)
      undef=-9.99e33
      do i=1,n*12
        sx(i)=undef
      enddo
      if(mon0.le.mon1)then
        nk=mon1-mon0+1
        nn=nk*n
        num=0
        do i=1,n
        do j=mon0,mon1
          num=num+1
          sx(num)=x(i,j)
        enddo
        enddo
      else
        nk=12+mon1-mon0+1
        nn=nk*n
        num=0
        do i=1,n
          do j=1,mon1
            num=num+1
            sx(num)=x(i,j)
          enddo
          do j=mon0,12
            num=num+1
            sx(num)=x(i,j)
          enddo
        enddo
      endif
      return
      end