c-----*----------------------------------------------------6---------7--
c     For a time vector series f(n,2), computing:
c       (1) mean vector in the high index years of x(n)
c       (2) mean vector in the low index years of x(n)
c       (3) composite difference between the mean vector of f(n,2) in high 
c                     index years of x(n) and its climatology average
c       (4) composite difference between the mean vector of f(n,2) in low 
c                     index years of x(n) and its climate average
c       (5) composite difference between the mean vectors of f(n,2) in high 
c                     and low index years of x(n)
c     input: n,x(n),f(n,2),coefh,coefl
c       n: the length of time series
c       x: control variable (index)
c       f: given series
c       coefh: control parameter for high index (i.e., high index years are 
c              those in which x(i) > coefh)
c       coefl: control parameter for low index (i.e., low index years are 
c              those in which x(i) < coefl)
c     output: fh(2),fl(2),dh(2),dl(2),dhl(2),tn(5)
c       fh: the mean vector of f in high index years of x(n)
c       fl: the mean vector of f in low index years of x(n)
c       dh: composite difference between the mean vector of f in high index years of x(n) 
c            and its climate mean (i.e., high index years minus cliamte mean).
c       dl: composite difference between the mean vector of f in low index years of x(n) 
c            and its climate mean (i.e., low index years minus cliamte mean). 
c       dhl: composite difference between the mean vectors of f in high and low index years
c             of x(n) (i.e., high minus low index years) 
c       tn(i,j): tn only equals 1. or 0. corresponding to significant difference or not. 
c           tn=1. indicates that the difference is significant.
c           tn=0. indicates that the difference is not significant.
c           tn(1,j)~tn(5,j) are corresponding to the 90%,95%,98%,99% and 99.9% confident levels. 
c           j=1: siginificant test for dh
c           j=2: siginificant test for dl
c           j=3: siginificant test for dhl
c     March 10, 2004 by Jianping Li.
      subroutine differhl1V(n,x,f,coefh,coefl,fh,fl,dh,dl,dhl,tn)
      dimension x(n),f(n,2)
      dimension fh(2),fl(2),dh(2),dl(2),dhl(2),tn(5,3)
      real avef(2),hn(n,2),ln(n,2) !work array
      do j=1,2
        call meanvar(n,f(1,j),avef(j),sf,vf)
      enddo
      do j=1,2
        nh=0
        nl=0
        do k=1,n
          if(x(k).gt.coefh)then
	    nh=nh+1
   	    hn(nh,j)=f(k,j)
	  endif
          if(x(k).lt.coefl)then
	    nl=nl+1
	    ln(nl,j)=f(k,j)
	  endif
        enddo
        call meanvar1(nh,hn(1,j),fh(j))
        call meanvar1(nl,ln(1,j),fl(j))
        dh(j)=fh(j)-avef(j)
        dl(j)=fl(j)-avef(j)
        dhl(j)=fh(j)-fl(j)
      enddo
      call diff_t_testV1(n,nh,f(1,1),f(1,2),hn(1,1),hn(1,2),tn(1,1))
      call diff_t_testV1(n,nl,f(1,1),f(1,2),ln(1,1),ln(1,2),tn(1,2))
      call diff_t_testV1(nh,nl,hn(1,1),hn(1,2),ln(1,1),ln(1,2),tn(1,3))
      return
      end
c-----*----------------------------------------------------6---------7--
c     Student's t-test for the difference between two vector means 
c     input: n,m,x(n,2),y(m,2),nt
c       n: number of x
c       m: number of y
c       x(n,2): vector series
c       y(m,2): vector series
c       nt: significant level, nt must be the one of 90,95,98,99 and 999
c     output: tn(5) 
c       tn: tn only equals 1. or 0. corresponding to significant difference or not. 
c           tn(1)~tn(5) are corresponding to the 90%,95%,98%,99% and 99.9% confident levels.
c     By Dr. Jianping Li, March 05, 2004 
      subroutine diff_t_testV1(n,m,x1,x2,y1,y2,tn)
      parameter(nn=10000)
      real x(n,2),y(m,2),tn(5),x1(n),x2(n),y1(m),y2(m)
      real ft(nn,5),ax(2),ay(2),wx(n,2),wy(m,2),wx1(n),wy1(m) !work array
      do i=1,n
        x(i,1)=x1(i)
        x(i,2)=x2(i)
      enddo
      do i=1,m
        y(i,1)=y1(i)
        y(i,2)=y2(i)
      enddo
      do i=1,5
        tn(i)=0.
      enddo
      do j=1,2
        call meanvar1(n,x(1,j),ax(j))
        call meanvar1(m,y(1,j),ay(j))
      enddo
      do j=1,2
        do i=1,n
          wx(i,j)=x(i,j)-ax(j)
        enddo
        do i=1,m
          wy(i,j)=y(i,j)-ay(j)
        enddo
      enddo
      do i=1,n
        wx1(i)=sqrt(wx(i,1)**2+wx(i,2)**2)
      enddo
      do i=1,m
        wy1(i)=sqrt(wy(i,1)**2+wy(i,2)**2)
      enddo
      call meanvar(n,wx1,ax1,sx1,vx1)
      call meanvar(m,wy1,ay1,sy1,vy1)
      sxy=(n*vx1+m*vy1)/(n+m-2.)
      sxy=sxy*(1./n+1./m)
      sxy=sqrt(sxy)
      if(sxy.eq.0.)return
      d1=ax(1)-ay(1)
      d2=ax(2)-ay(2)
      dxy=sqrt(d1**2+d2**2)
      ta=dxy/sxy
      nm=n+m-2
      call t_table(ft(1,1),ft(1,2),ft(1,3),ft(1,4),ft(1,5))
      do i=1,5
        if(ta.ge.ft(nm,i))then
          tn(i)=1.
        endif
      enddo
      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
c-----*----------------------------------------------------6---------7--
c     Computing the mean ax 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 and sx
c       ax: the mean value of x(n)
c     By Dr. LI Jianping, May 6, 1999.
      subroutine meanvar1(n,x,ax)
      dimension x(n)
      ax=0.
      vx=0.
      sx=0.
      do 10 i=1,n
        ax=ax+x(i)
  10  continue
      ax=ax/float(n)
      return
      end
c-----*----------------------------------------------------6---------7--
c     t-distribution, i.e., student's distribution
c     t table with two-tailed (right and left tails) probabilities
c       P(|t|>=ta)=a
c       where a (alpha) significance level and (1-a)*100% is confindence level.
c     t90: t-distribution test at 90% confidence level.
c     t95: t-distribution test at 95% confidence level.
c     t98: t-distribution test at 98% confidence level.
c     t99: t-distribution test at 99% confidence level.
c     t999: t-distribution test at 99.9% confidence level.
c     By Dr. Jianping Li, January 5, 2000.
      subroutine t_table(ft90,ft95,ft98,ft99,ft999)
      parameter(n=10000)
      dimension ft90(n),ft95(n),ft98(n),ft99(n),ft999(n)
      dimension n90(30),n95(30),n98(30),n99(30),n999(30)
      data n90/ 6314,2920,2353,2132,2015,1943,1895,1860,1833,1812,
     *          1796,1782,1771,1761,1753,1746,1740,1734,1729,1725,
     *          1721,1717,1714,1711,1708,1706,1703,1701,1699,1697/
      data n95/12706,4303,3182,2776,2571,2447,2365,2306,2262,2228,
     *          2201,2179,2160,2145,2131,2120,2110,2101,2093,2086,
     *          2080,2074,2069,2064,2060,2056,2052,2048,2045,2042/
      data n98/31821,6965,4541,3747,3365,3143,2998,2896,2821,2764,
     *          2718,2681,2650,2624,2602,2583,2567,2552,2539,2528,
     *          2518,2508,2500,2492,2485,2479,2473,2467,2462,2457/
      data n99/63657,9925,5841,4604,4032,3707,3499,3355,3250,3169,
     *          3106,3055,3012,2977,2947,2921,2898,2878,2861,2845,
     *          2831,2819,2807,2797,2787,2779,2771,2763,2756,2750/
      data n999/636619,31598,12941,8610,6859,5959,5405,5041,4781,
     *     4587,4437,4318,4221,4140,4073,4015,3965,3922,3883,3850,
     *          3819,3792,3767,3745,3725,3707,3690,3674,3659,3646/
      do 10 i=1,30
        ft90(i)=n90(i)/1000.
        ft95(i)=n95(i)/1000.
        ft98(i)=n98(i)/1000.
        ft99(i)=n99(i)/1000.
        ft999(i)=n999(i)/1000.
  10  continue
      do 20 i=31,40
        fi=(i-30.)/10.
        ft90(i)=1.697-(1.697-1.684)*fi
        ft95(i)=2.042-(2.042-2.021)*fi
        ft98(i)=2.457-(2.457-2.423)*fi
        ft99(i)=2.750-(2.750-2.704)*fi
        ft999(i)=3.646-(3.646-3.551)*fi
  20  continue
      do 30 i=41,60
        fi=(i-40.)/20.
        ft90(i)=1.684-(1.684-1.671)*fi
        ft95(i)=2.021-(2.021-2.000)*fi
        ft98(i)=2.423-(2.423-2.390)*fi
        ft99(i)=2.704-(2.704-2.660)*fi
        ft999(i)=3.551-(3.551-3.460)*fi
  30  continue
      do 40 i=61,120
        fi=(i-60.)/60.
        ft90(i)=1.671-(1.671-1.658)*fi
        ft95(i)=2.000-(2.000-1.980)*fi
        ft98(i)=2.390-(2.390-2.358)*fi
        ft99(i)=2.660-(2.660-2.617)*fi
        ft999(i)=3.460-(3.460-3.373)*fi
  40  continue
      do 50 i=121,n
        fi=(i-120.)/(n-120.)
        ft90(i)=1.658-(1.658-1.645)*fi
        ft95(i)=1.980-(1.980-1.960)*fi
        ft98(i)=2.358-(2.358-2.326)*fi
        ft99(i)=2.617-(2.617-2.576)*fi
        ft999(i)=3.373-(3.373-2.291)*fi
  50  continue
      return
      end