c-----*----------------------------------------------------6---------7--
c     For a time series f(n), computing:
c       (1) mean in the nc strongest index years of x(n)
c       (2) mean in the nc weakest index years of x(n)
c       (3) composite difference between the mean of f(n) in nc strongest 
c                     index years of x(n) and its climatology average
c       (4) composite difference between the mean of f(n) in nc weakest 
c                     index years of x(n) and its climate average
c       (5) composite difference between the means of f(n) in nc strongest 
c                     and weakest index years of x(n)
c     input: n,x(n),f(n),nc
c       n: the length of time series
c       x: control variable (index)
c       f: given series
c       nc: number of the strongest positive (or negative) values of x
c     output: fh,fl,dh,dl,dhl,tn(5)
c       fh: the mean of f in nc strongest index years of x(n)
c       fl: the mean of f in nc weakest index years of x(n)
c       dh: composite difference between the mean of f in nc strongest index years of x(n) 
c            and its climate mean (i.e., strong index years minus cliamte mean).
c       dl: composite difference between the mean of f in nc weakest index years of x(n) 
c            and its climate mean (i.e., weak index years minus cliamte mean). 
c       dhl: composite difference between the means of f in nc strongest and nc weakest index years
c             of x(n) (i.e., strong minus weak index years) 
c       tn(i,j): tn only equals 2., -2., 1., -1. or 0. corresponding to significant difference or not.
c           tn=2. indicates that the difference is positive and significant.
c           tn=-2. indicates that the difference is negative and significant.
c           tn=1. indicates that the difference is positive but not significant.
c           tn=-1. indicates that the difference is negative but not significant.
c           tn=0. indicates the difference is zero.
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 differencehl2(n,x,f,nc,fh,fl,dh,dl,dhl,tn)
      dimension x(n),f(n),tn(5,3),ip(n)
      real hn(nc),ln(nc) !work array
      call meanvar(n,f,avef,sf,vf)
      fh=0.
      fl=0.
      call labelsort(n,1,1,x,ip)
      do k=1,nc
        hn(k)=f(ip(k))
        ln(k)=f(ip(n-k+1))
      enddo
      call meanvar1(nc,hn,fh)
      call meanvar1(nc,ln,fl)
      dh=fh-avef
      dl=fl-avef
      dhl=fh-fl
      call diff_t_test(nc,n,hn,f,tn(1,1))
      call diff_t_test(nc,n,ln,f,tn(1,2))
      call diff_t_test(nc,nc,hn,ln,tn(1,3))
      return
      end
c-----*----------------------------------------------------6---------7--
c     Lable sort method for a two-dimensional series
c     input: n,m,nl,x(n,m)
c       n: number of time series
c       m: number of elements of x 
c       nl: specified colucmn as label (1<=nl<=m)
c       x(n,m): the raw series
c     output: ip(n)
c       ip(n): label array, indicating the orignal position of the order series in x(n)
c          the order series is arrayed from big to small
c     March 09, 2004 by Jianping Li.
      subroutine labelsort(n,m,nl,x,ip)
      dimension x(n,m),ip(n)
      do i=1,n
        ip(i)=i
      enddo
      n1=n-1
      do 30 i=1,n1
        max=i
        jj=i+1
        do j=jj,n
          if(x(ip(j),nl).gt.x(ip(max),nl))then
            max=j
          endif
        enddo
        itm=ip(i)
        ip(i)=ip(max)
        ip(max)=itm
  30  continue
      return
      end
c-----*----------------------------------------------------6---------7--
c     Student's t-test for the difference between two means 
c     input: n,m,x(n),y(m)
c       n: number of x
c       m: number of y
c       x(n): series
c       y(m): series
c     output: tn(5) 
c       tn: tn only equals 2, -2., 1, -1. or 0. corresponding to significant difference or not. 
c           tn=2. indicates that the difference is positive and significant.
c           tn=-2. indicates that the difference is negative and significant.
c           tn=1. indicates that the difference is positive but not significant.
c           tn=-1. indicates that the difference is negative but not significant.
c           tn=0. indicates the difference is zero.
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_test(n,m,x,y,tn)
      parameter(nn=10000)
      real x(n),y(m),tn(5)
      real ft(nn,5) !work array
      do i=1,5
        tn(i)=0.
      enddo
      call meanvar(n,x,ax,sx,vx)
      call meanvar(m,y,ay,sy,vy)
      sxy=(n*vx+m*vy)/(n+m-2.)
      sxy=sxy*(1./n+1./m)
      sxy=sqrt(sxy)
      if(sxy.eq.0.)return
      dxy=ax-ay
      if(dxy.gt.0.)then
        sn=2.
        do i=1,5
          tn(i)=1.
        enddo
      else
        sn=-2.
        do i=1,5
          tn(i)=-1.
        enddo
      endif
      ta=abs(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)=sn
        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
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