      PROGRAM ep
C
C PGM TO PERFORM HOMOGENEITY TESTS ON climate time series using 
C LINEAR REGRESSION with unspecified changepoints.
C Potential changepoint is identified as the
c time series place with minimum RSS in Fratio, or with student's t, 
c at 0.05
c then verified as significant using MRPP test at .05 probability.
C with window. Final offset is calculated with a window +-10.
c
c Bug fixed 9/6/95
c
C D. EASTERLING - T. Peterson  12/8/93
C   Last update/modification   9/7/95
C MAIN PROGRAM
      real*8 A(500),B(500),DIFF1(500),falpha2,talpha
      real*8 toma(500),tomb(500),tomd(500),tomdtcorr
      real*8 MISS,Minsumsq0,minsumsq1,df1a,xtemp(500),ytemp(500)
      real*8 minsumsq,n,df0,df1,falpha,ta,tb
      real*8 a1,b1,a2,b2,a3,b3,t,corr2
      real*8 savef(50),f,talpha2,atemp(500),numer,denom	
      INTEGER I,J,NUM,L,Q,chgpt,flags(500),change(50),length,yr(500)
      integer beg,end,flag,fchange(50),fchange2(50,2),iq2,endflag
      integer tempyr(500),igoodnuf,choice
      integer tyr,tid,id(500,2),lagid,tempflag(500),chgpt2
      integer id2,lagid2
c      character*6 tid,id(500),lagid
      CHARACTER*50 NAMe1,NAME2,name3
      WRITE(6,125)
 125  FORMAT(1X,'enter file name for homogeneity testing')
      READ(5,135)NAME1
 135  FORMAT(A50)
      WRITE(6,145)

c	we are analyzing tmax or tmin

	write (*,*) ' '
	write (*,*) 'Please Select '
	write (*,*) ' '
	write (*,*) '1. Maximum Temperature'
	write (*,*) '2. Minimum Temperature' 
	write (*,*) ' '
	read (*,*) choice
	write (*,*) ' '
	
   

 145  FORMAT(1X,'Output Files?, 1st for stats, 2nd for adj series')
      READ(5,135)NAME2
      read(5,135)name3
       write(6,*)'Enter 1 if data are degrees F and are to be converted'
       write(6,*)'   or 2 if they are not to be converted'
      read(5,*) iconvert
      OPEN(20,FILE=name1)
      OPEN(42,FILE=name2)
      open(52,file=name3)
	if (choice .eq. 1) then
      	  open(62,file='tmax.error.file.HII.F')
	else
	  open(62,file='tmin.error.file.HII.F')
	endif
C
C MAIN PROGRAM LOOP, ONE TIME THROUGH FOR EACH TIME SERIES
C
C
      talpha=10.
      falpha=10.
      FLAG=1
      endflag=0
      NUM=0
      lagid=999999
      MISS=-99.990
      L=0
c
       goto 5
 4    num=1
      lagid=tid
      lagid2=id2
      id(num,1)=tid
      id(num,2)=id2
      yr(num)=tyr
      a(num)=ta
      b(num)=tb
      toma(num)=ta
      tomb(num)=tb
      tomd(num)=miss
      if(toma(num).gt.-700.and.tomb(num).gt.-700)
     +tomd(num)=ta-tb

 5     num=num+1
C
C READ TWO TIME SERIES 
C
c various read statements for differenct data sets, choose one
c a is the candidate series array, b is the reference series array
c  input file generally looks like:  station id number, year, 
c  candidate station data series, reference series value, difference
c  (difference or ratio is calculated internally so it is not necessary as
c  input)
c        READ(20,197,end=20) tcountry,tid,tyr,ta,tb,canddt,
c       read(20,*,end=20)tid,tyr,ijunk,ta,tb,rjunk  
       read(20,*,end=20)tid,id2,tyr,ijunk,ta,tb,rjunk  
c       read(20,198,end=20) tid,tyr,ta,tb
c 197    format(i3,1x,i7,1x,i4,4(f10.3),2(f7.3))
  198    format(i7,1x,i4,2(f10.3))
c
c convert from deg f to deg C
c
       if (iconvert.eq.1) then
        if (ta.gt.miss+1) then
        ta=(ta-32.)*.5556
        end if
        if (tb.gt.miss+1) then
        tb=(tb-32.)*.5556
        end if
       end if
      if ((tid.ne.lagid.or.id2.ne.lagid2)
     +.and.(lagid.ne.999999)) goto 30
      tomdtcorr=corr2
      lagid=tid
      lagid2=id2
      id(num,1)=tid
      id(num,2)=id2
      yr(num)=tyr
      a(num)=ta
      b(num)=tb
      toma(num)=ta
      tomb(num)=tb
      tomd(num)=miss
      if(toma(num).gt.-700.and.tomb(num).gt.-700)
     +tomd(num)=ta-tb
      goto 5
 20   endflag=1
 30   num=num-1
      L=L+1
c      print *,'l and num',l,num
      ifirst=0
      ipass=0
      do j=1,num
       flags(j)=0
      end do
c set temporary array for testing of cand series
        do i1=1,num
         atemp(i1)=a(i1)
        end do

c  check to see if there are enough years of reference series for the
c  station and if the correlation is good enough
         call good_nuf(miss,a,b,igoodnuf,num)
c    this subroutine may shorten the reference series to shorten if early
c    years weren't very highly correlated

	if(igoodnuf.eq.1)then !not good enough station data to use
c  set adjusted to missing, then go to the write adjusted file out spot
c   from there it will loop back to next read
	  do i11=1,num
		a(i11)=miss 
               	write(52,360)id(i11,1),id(i11,2)
     +,yr(i11),a(i11),toma(i11)
 360     format(1x,i8,i3,1x,i4,1x,2f8.2)


	  enddo
	  go to 98
	endif


c
c form difference series
c
 50   CALL DIFFER(ATEMP,B,NUM,MISS,FLAG,DIFF1)
c 
c reset flags array, if a discont is found, set =4
c to check in changepoints1 subroutine.  Flags array keeps track
c of where the discontinuities are.
c 
      talpha=10.
      falpha=10.
      do j=1,num
        if (flags(j).eq.2) then
        flags(j)=4
        else if (flags(j).ne.4) then 
        flags(j)=0
        end if
      end do
c   The following lines set begin and end  points of diff series
c   to be tested
c    
       q=0
  54   i=0
       beg=0
       end=0
  55   i=i+1
   	if ((flags(i).eq.0).and.(beg.eq.0)) then
c  set begin point
         beg=i
         goto 56 
        end if
        if ((beg.eq.0).and.(i.eq.num-3)) goto 95
        goto 55
 56      i=i+1
        if (((flags(i).eq.2).or.((flags(i).eq.3).and.(ipass.gt.1)))
     c  .or.((flags(i).eq.1)
     c     .and.(end.eq.0))) then
c  set end point
         end=i-1
         goto 57
         end if
        if ((i.eq.num).and.(end.eq.0)) then
c set end for first time through
         end=num
         goto 57
        end if
        goto 56
 57      j=0
       if (end-beg.lt.11) then
c 
c time series is too short to test
c
        t=0.0
        f=0.0
        goto 93
       end if
c
c set test arrays, ytemp is difference series, xtemp is time 
c                                              (e.g. yrs 1950-1990)
c
       do i=beg,end
        if (diff1(i).gt.miss+1) then
        j=j+1
        ytemp(j)=diff1(i)
        xtemp(j)=float(i)
        tempflag(j)=flags(i)
        tempyr(j)=yr(i)
        if (tempflag(j).eq.4) then
c         print *,beg,end,tempyr(j),ipass
        end if
        end if
       enddo
        if (j.lt.11) then
c this section of time series is too small to test, goto 93
          t=0.0
          f=0.0
         goto 93
        end if
 175        length=j
c
c set critical values for t and f tests
c
        call alpha(length,falpha,falpha2,talpha,talpha2)
c
c call subroutines for 0 and 1 possible changepoint
c

       CALL chngpts0(ytemp,xtemp,j,MISS,minsumsq0,a1,b1)
       minsumsq=minsumsq0
       call chngpts1(ytemp,xtemp,j,miss,chgpt,minsumsq,minsumsq1,
     c a2,b2,a3,b3,tempflag,l,chgpt2,tempyr,ipass)
c
c
       if (minsumsq1.eq.minsumsq0) then
c if so then the changpoints1 subroutine didn't execute the 
c regression since other potential changpoints where too close
c and ytemp series was too small to get past the +-4 places 
c requirement for how close together potential disconts can be.
         f=0.0
         t=0.0
         goto 93
        end if
       minsumsq=minsumsq1
       n=float(j)
       df0=n-5.
       df1=n-8.
       df1a=n-4.
c       u=((minsumsq0-minsumsq1)/3.)/(minsumsq1/(df1a))
c
c u likelihood ratio is from Solow (1987) JCAM article on 2-phase
c regression, it has an F distribution, not used !!!
c
        if (df1.gt.0.0) denom=minsumsq1/df1
        if (df0-df1.gt.0.0) numer=(minsumsq0-minsumsq1)/(df0-df1)
       if (denom.gt.0.0) then
	f=numer/denom
	else
	  f=0.0
       end if
c
c f is the f-ratio from Lucie Vincent's method (Canadian Climate Centre)
c call ttest to test for sig difference in means
c t-test makes sure that a trend in the difference series is picked
c up as a series of discontinuities
c
        call ttest(ytemp,j,chgpt2,t)
c
c write to test file (52) prior to checking for significance
c also set potential changepoint = 3, this will change to 2 if found
c to be significant
c
c        write(52,300)l,chgpt,u,f,tt,ttyhat,beg,end
	flags(chgpt)=3
	 if ((t.ge.talpha).or.(f.ge.falpha)) then
c then we have a winner, set changepoint from potential (=3) to 
c confirmed (=2) and ifirst=1 and try again
c  
c test print for debugging
c
c        if (l.eq.6) print *,chgpt, ' chgpt ', beg,end
c
         flags(chgpt)=2
        q=q+1
        change(q)=chgpt
        savef(q)=prob1
      else
         t=0.0
         f=0.0
        end if
  93	if ((t.lt.talpha).and.(f.lt.falpha)) then
c this part of series is homogeneous or too small, set all flag (i)
c between begin and end to 1, indicating it has been tested
        do i=beg,end
         if (flags(i).eq.0) flags(i)=1
        enddo
       if (end.eq.num) goto 95
       end if
       goto 54
 95    if (q.gt.0) then
c
c create temp cand series, call window to verify, then test adjusted
c temp series again
c
	ipass=ipass+1
	call window(atemp,diff1,change,savef,q,num,fchange,iq,l,miss)
        if ((iq.eq.0).and.(ipass.eq.1)) then
          do i=1,num
          if (flags(i).eq.1) flags(i)=0
          enddo
          goto 54
        else if (iq.gt.0) then
c keep verified disconts in fchange2 array for passing to final adjust
         do i9=1,iq
c         print *,fchange(i9),ipass
           fchange2(i9+iq2,1)=fchange(i9)
           fchange2(i9+iq2,2)=ipass
           ifirst=2
         end do
         iq2=iq2+iq
        else if ((iq.eq.0).and.(ipass.gt.1)) then 
           goto 91         
        end if
c
c temp cand series has just been adjusted using windows so goto 50
c and check for more disconts in adjusted time series, if ipass2 is
c 0 then no further disconts were found, so write out adjusted series
c and go to next time series
c
c	print *,ipass
c	write(42,*)ipass
        goto 50
	end if
c
c write adjusted candidate series to file 52 
c first re-create difference series to pass to adjust
c	
  91      CALL DIFFER(A,B,NUM,MISS,FLAG,DIFF1)
  	call adjust(a,b,diff1,fchange2,num,iq2,yr,l,miss,id,
     +toma,tomb,tomd,tomdtcorr)
c reset iq2 to zero
	iq2=0
C  
c      if (l.gt.20) goto 99  
  98	continue  !loops here from writing bad series post good_nuf
       if (endflag.eq.1) goto 99        
      goto 4     
 99   continue
      write(*,*)'homogeneity test complete'

      END
c
c *************************
c end of main program
c *************************
c        
	subroutine window(a,diff1,change,savef,q,num,fchange,iq,l
     c   ,miss)
c subroutine to verify discontinuities using multiresponse
c permutation procedures (mrpp) see Mielke, et al. (1981) MWR pp. 120-126
c subroutine also applies adjustments to candidate series during testing 
c uses difference series means for adjustment factor
c
	real*8 a(500),diff1(500),savef(50),miss
	real*8 dy(500),prob,ptest,lprob
	real*8 ftemp,x1,mean1,x2,mean2,offset
	integer q,change(50),beg1,end1,beg2,end2,chgpt
        integer beg,end,fchange(50),iq
        q=q+2
	ptest=0.05
       savef(q)=0.0
       savef(q-1)=0.0
       change(q)=1
       change(q-1)=num
c sort changepoints to begin adjusting at end of time series
      do i1=1,q-1
       do i2=i1+1,q
	if (change(i1).lt.change(i2)) then
          itemp=change(i1)
          change(i1)=change(i2)
          change(i2)=itemp
          ftemp=savef(i1)
          savef(i1)=savef(i2)
          savef(i2)=ftemp
     	 end if
	end do
c end do i2
      enddo
c end do i1 
	iq=0
c 
c end of sort and begin testing changepoints
c
	do i1=1,q-2
          beg1=change(i1+1)+1
	  end1=change(i1)
          beg2=change(i1+2)
          end2=change(i1+1)
         
          i5=0
          lprob=1.0
          beg=change(i1+1)-11
          end=change(i1+1)+12
          if (beg.lt.beg2) beg=beg2
          if (end.gt.end1) end=end1
          i5=0
          do i4=beg,end
           if (diff1(i4).gt.miss+1) then
             i5=i5+1
             if (i4.eq.beg1) chgpt=i5
             dy(i5)=diff1(i4)
           end if
         end do
c make sure each side of time series before and after chgpt has 2
c or more values, otherwise skip out and set prob to 1.0 
         if ((chgpt.ge.4).and.(i5-chgpt.ge.3)) then
c end do i4
c call multiresponse permutation procedure subroutine
c
 	call mrppt(chgpt,i5,dy,prob)
c        print *,'from mrpp ',prob,beg1,l
        else
         prob=1.0
        end if
        if (prob.le.ptest) then
c adjust series with window concept if prob le ptest
        mean1=0.0
          x1=0.0
         do i2=beg,change(i1+1)
          if (diff1(i2).gt.miss+1) then
	   mean1=mean1+diff1(i2)
	   x1=x1+1.0
	  end if
	 end do
c end do i2
          mean2=0.0
          x2=0.0
 	 do i2=change(i1+1)+1,end
	  if (diff1(i2).gt.miss+1) then
           mean2=mean2+diff1(i2)
           x2=x2+1.0
          end if
         end do 
c end do i2
         if ((x1.gt.0.0).and.(x2.gt.0.0)) then
          mean1=mean1/x1
         mean2=mean2/x2
         offset=mean2-mean1
c         print *,mean1,mean2,offset,end2,change(i1+1) 
         iq=iq+1
          fchange(iq)=change(i1+1)  
	 do i2=1,end2
          if (A(i2).gt.miss+1) then
             a(i2)=A(i2)+offset
          end if
         end do
c end do i2
         end if
         end if
        end do
c end do q
        return
        end
c
c  final adjustment subroutine
c
       subroutine adjust(a,b,diff1,change,num,q,yr,l,miss,id,
     +toma,tomb,tomd,tomdtcorr)
	real*8 a(500),diff1(500),x1,mean1,x2,mean2,offset,b(500)
        real*8 miss,toma(500),tomb(500),tomd(500),tomdtcorr,cor
        real*8 adt(500),bdt(500),atmp(500),varx,vary,n,avgx,avgy
        real*8 sumx,sumy
	integer q,change(50,2),l,beg1,end1,beg2,end2,beg,end
        integer id(500,2),change2(50,2),yr(500)
c	print *,'adjust q=',q
        if (q.eq.0) goto 99 !       
c
c create temp and final a() time series
c
        do i1=1,num
         atmp(i1)=a(i1)
        end do
c if q=0 then do not adjust but write out anyway
        q=q+2
       change(q,1)=1
       change(q-1,1)=num
       change(q,2)=100
       change(q-1,2)=100
c sort changepoints to begin adjusting at end of time series
      do i1=1,q-1
       do i2=i1+1,q
	if (change(i1,1).lt.change(i2,1)) then
          itemp1=change(i1,1)
          itemp2=change(i1,2)
          change(i1,1)=change(i2,1)
          change(i1,2)=change(i2,2)
          change(i2,2)=itemp2
          change(i2,1)=itemp1
    	 end if
	end do
c enddo i2
      enddo 
c end do i1
c
c filter out changepoints closer than 5 places apart
c
         iq=1
c set last year of time series (e.g. 100)
         change2(1,1)=change(1,1)
         change2(1,2)=100
         i1=1
 5       i1=i1+1
         
         if (change(i1,1).lt.change(i1+1,1)+5) then
c write to error file
         write(62,*)'Error, discont too close at
     c ',id(1,1),id(1,2),yr(change(i1,1)),yr(change(i1+1,1))
c if too close determine which was found first 
	   if (change(i1,2).le.change(i1+1,2)) then
              iq=iq+1
              change2(iq,1)=change(i1,1)
              change2(iq,2)=change(i1,2)
              i1=i1+1
           else if (change(i1,2).ge.change(i1+1,2)) then
               iq=iq+1
               change2(iq,1)=change(i1+1,1)
	       change2(iq,2)=change(i1+1,2)
               i1=i1+1
           end if
         else
c if far enough apart then copy current change(i1,1) to change2
                iq=iq+1
                change2(iq,1)=change(i1,1)
                change2(iq,2)=change(i1,2)
         end if
        if (i1.lt.q-1) goto 5
	iq=iq+1
c set first year of time series (e.g. 1)
         change2(iq,1)=change(q,1)
         change2(iq,2)=1

	do i1=1,iq-2
	  end1=change2(i1,1)
          beg1=change2(i1+1,1)+1
	  end2=change2(i1+1,1)
          beg2=change2(i1+2,1)
c
c create dt/dT series for cand and ref subsections from previous
c discont to next discont, excluding year of current discont, then
c call correlation
c see Peterson and Easterling (1994) Int. J. Clim. pp. 671-679 
c for discussion of dt/dT series.
c
          i3=0
          do ii2=beg2,end1-1
           if (ii2.eq.end2)atmp(ii2)=miss
           if ((atmp(ii2+1).gt.miss+1).and.(atmp(ii2).gt.miss+1)
     c .and.(b(ii2+1).gt.miss+1).and.(b(ii2).gt.miss+1))then
            i3=i3+1
            adt(i3)=a(ii2)-a(ii2+1)
            bdt(i3)=b(ii2)-b(ii2+1)
           end if
          end do 
         atmp(end2)=a(end2)
      call CORREL(adt,bdt,i3,miss,COR,VARX,VARY,N,AVGX,AVGY,SUMX,SUMY)
c
c if r is less than 0.8 then skip out, not strong enough relationship
c between the subsections of cand and ref series.
c  somewhere along this gradiant, as the correlation drops between the
c  ref and candidate, spurious things will occur in the reference series.
c  at what point is that reached.  Tom Peterson tested this in the summer
c  of 1994 and determined that it was at an r (dT/dt series) of 0.79.  So
c  we will round up to an even 0.80.
c
c These time series must now be filtered out of the input data set
c
c      if (cor.lt.0.8) goto 97     
c
c calculate offset as diff of means of diff series before and
c after changepoint, with other changepoints as boundaries of window
c
	  mean1=0.0
          x1=0.0
	 beg=end2-11
         if (beg.lt.beg2) beg=beg2
         do i2=beg,end2
          if (diff1(i2).gt.miss+1) then
	  mean1=mean1+diff1(i2)
	  x1=x1+1.0
	  end if
	 end do
c end do i2
          mean2=0.0
          x2=0.0
	 end=beg1+12
         if (end.gt.end1) end=end1
 	 do i2=beg1,end
	  if (diff1(i2).gt.miss+1) then
           mean2=mean2+diff1(i2)
           x2=x2+1.0
          end if
         end do 
c end do i2
         if ((x1.gt.0.0).and.(x2.gt.0.0)) then
         mean1=mean1/x1
         mean2=mean2/x2
         offset=mean2-mean1
c 
c if discontinuity is larger than +- 3.0C for annual
c or +- 6.0C for seasonal then delete 
c time series from final data set
c
         if (abs(offset).gt.6.0) goto 999
c
c adjust series
c
	 do i2=1,end2
          if (A(i2).gt.miss+1) then
             a(i2)=A(i2)+offset
          end if
         end do
c end do i2
         end if
c
c write statistics and changepoints, etc to file 42
c
c        print *,change2(i1+1,1),change2(i1+1,2)
c for debugging:

c for debugging:
cccc

c for debugging:
cc        write(4300)l,id(1),change2(i1+1,1),change2(i1+1,2),
c     c  yr(change2(i1+1,1)),offset,
c     c  yr(beg2),yr(end1),tomdtcorr,beg2,end1,cor
c 300   format (1x,i4,1x,i8,1x,i4,1x,i2,1x,i4,1x,f8.2,1x,i4,1x,i4,
c     +f8.2, 1x,i4,1x,i4,1x,f4.2)

       write(42,300)l,id(1,1),id(1,2),
     +yr(change2(i1+1,1)),offset,yr(beg2),yr(end1),cor
 300   format (1x,i4,1x,i8,i3,1x,i4,1x,f8.2,1x,i4,1x,i4,f6.2)


cc
cc  the output file for
cc
cc Example output from stats file, each line is one discontinuity
c
c  L        id    yr     offset  byr  eyr    dtcor
c 3400     3520  1956     0.40   1946 1992    0.88
c
c where: L     loop counter, Lth station processed
c       id    station id
c       yr    the actual year the discontinuity was found
c      offset the adjustment factor, which is the difference between the means
c             of the difference series before and after the discontinuity
c       byr   begin year of the subsection of time series tested
c       eyr   end year of the subsection of the time series tested
c       dtcor is the correlation between 1st difference series for that section
c
 97    end do
c end do iq
c
c write final adjusted series to unit 52
c
c   write out heading at top of columns to keep them straight
c	write(52,*)'                   adjusted   raw     ref',
c     +'  origdif finaldif adjustment dtcorr'
c
c  If the early part of the reference and candidate series is not
c correlated enough, (from sub good_nuf) but becomes so later on, 
c make the early part of the adjusted candidate series missing.
c
	do i11=1,num
	  if(b(i11).lt.miss+1)then !if no early reference series
	   a(i11)=miss
	  else
	    go to 99
	  endif
	enddo


 99     do i11=1,num
	tomg=-9999
	tomh=-9999
	if(a(i11).gt.-700.and.tomb(i11).gt.-700)tomg=a(i11)-tomb(i11)
	if(a(i11).gt.-700.and.toma(i11).gt.-700)tomh=a(i11)-toma(i11)
c  hence tomh=adjusted minus raw=adjustment
c        tomg=adjusted-ref=final (adjusted) difference series
c  in debug mode:
c          write(52,360)l,id(i11),yr(i11),a(i11),toma(i11),tomb(i11),
c     +  tomd(i11),tomg,tomh,tomdtcorr
c 360     format(1x,i6,1x,i8,1x,i4,1x,7f8.2)
c note for future changes::::there is another write to 52 statement
c right after the call good_nuf
          write(52,360)id(i11,1),id(i11,2),yr(i11),a(i11),toma(i11)
 360     format(1x,i8,i3,1x,i4,1x,2f8.2)

c   standard output now is:
c	station id, year, adjusted station value, raw station value

c  from debug mode output:
c
cExample output from adjusted file
c
c L       id    ayr    adj      raw     ref   origdif  finaldif offset dtcor    
c3400     3520 1907    18.56   15.97   19.12   -3.14   -0.56    2.58    0.00
c3400     3520 1908    18.31   15.73   18.43   -2.70   -0.12    2.58    0.00
c3400     3520 1909    18.22   15.64   17.59   -1.95    0.63    2.58    0.00
c3400     3520 1910    19.09   16.51   18.29   -1.78    0.80    2.58    0.00
c
cwhere: L     loop counter, Lth station processed
c       id    station id
c     ayr   year in time series
c       adj   adjusted value
c       raw   raw (unadjusted) value
c       ref   reference series value
c       origdif  raw minus ref
c       finaldif  adj minus ref
c       offset  adj minus raw

        end do
 999	return
	end
c 
c routine to create difference series
c
      SUBROUTINE DIFFER(A,B,NUM,MISS,FLAG,C)
      real*8 A(500),B(500),C(500),MISS
      INTEGER NUM,FLAG
      IF (FLAG.EQ.1) THEN
C DO TEMPERATURE
      DO 10 I=1,NUM
       IF ((A(I).gt.MISS+1).AND.(B(I).gt.MISS+1)) THEN
         C(I)=A(I)-B(I)
       ELSE
         C(I)=MISS
       END IF
 10    continue
       END IF
      IF (FLAG.EQ.2) THEN
C DO PRECIP
      DO 20 I=1,NUM
       IF ((A(I).gt.MISS+1).AND.(B(I).gt.MISS+1)) THEN
         C(I)=A(I)/B(I)
       ELSE
         C(I)=MISS
       END IF
 20    continue
       END IF
       RETURN
      END
c
	subroutine chngpts0(y,x,n,miss,minsumsq0,a1,b1)
          real*8 x(500),y(500),miss,minsumsq0,a1,b1,sumsq
          integer n
          call regress(x,y,miss,n,sumsq,a1,b1)
          minsumsq0=sumsq
          return
          end
c
	subroutine chngpts1(b,a,n,miss,chngpt1,minsumsq,minsumsq1,
     c      a12,b12,a22,b22,tempflag,l,chngpt2,tempyr,ipass)
         integer beg1,end1,beg2,end2,chngpt1,n,tempyr(500)
         integer skip,tempflag(500),l,chngpt2
         real*8 a(500),b(500),x(500),y(500),minsumsq,rsumsq
         real*8 a1,b1,a2,b2,sumres1,sumres2,miss,minsumsq1
         real*8 a12,b12,a22,b22
	 minsumsq1=minsumsq
         do 10 i=5,n-5
          beg1=1
	  skip=0
          icz=i-5
          if(icz.le.0)icz=1
c 
c skip to end of loop if year being checked is too close to
c another discont.(+- 5).
c
          do i9=icz,i+6
           if (tempflag(i9).eq.4) then
c           print *,'tempyr',tempyr(i9),i9,ipass
           goto 10
           end if
          end do
          end1=i
          i1=0
          do 20 j=beg1,end1
           i1=i1+1
           x(i1)=a(j)
           y(i1)=b(j)
 20       continue
          call regress(x,y,miss,i1,sumres1,a1,b1)
          beg2=end1+1
          end2=n
           i2=0
          do 30 j=beg2,end2
           i2=i2+1
           x(i2)=a(j)
           y(i2)=b(j)
 30       continue
          call regress(x,y,miss,i2,sumres2,a2,b2)
          rsumsq=sumres1+sumres2
          if (rsumsq.lt.minsumsq1) then
           minsumsq1=rsumsq
           chngpt1=idint(a(end1))
           chngpt2=end1
           a12=a1
           b12=b1
           a22=a2
           b22=b2
          end if
 10       continue
          return
          end
c
c
	Subroutine Regress(x,y,miss,num,sumres,a1,b1)
c **2 linear regression routine
        real*8 x(500),y(500),miss,sumres,resid(500),yhat
        real*8 b1,b21,b11,a1,xbar,ybar,n1
        integer num,i
	xbar=0.0
        ybar=0.0
	b1=0.0
        b11=0.0
        b21=0.0
        n1=0.0
        a1=0.0
	do 10 i=1,num
         if ((y(i).gt.miss+1).and.(x(i).gt.miss+1)) then
         n1=n1+1.0
         ybar=ybar+y(i)
         xbar=xbar+x(i)
         end if
 10     continue
        if (n1.gt.0.0) then
        xbar=xbar/n1
        ybar=ybar/n1
        end if
	do 15 i=1,num
         if ((x(i).gt.miss+1).and.(y(i).gt.miss+1)) then
         b11=b11+((x(i)-xbar)*(y(i)-ybar))
         b21=b21+((x(i)-xbar)**2)
         end if
 15     continue
        if (b21.gt.0.0) then 
        b1=b11/b21
        else
        b1=0.0
        end if
        a1=ybar-b1*xbar
        sumres=0.0
c calculate residuals, and residual sum of squares
	do 20 i=1,num
         if ((x(i).gt.miss+1).and.(y(i).gt.miss+1)) then
           yhat=a1+b1*x(i)
          resid(i)=y(i)-yhat
	  sumres=sumres+resid(i)*resid(i)
         end if
 20     continue
        return
        end
        SUBROUTINE TTEST(A,N,end1,TT)
C ROUTINE TO CALCULATE STUDENT'S T TEST
C FOR USE IN HOMOGENEITY ANALYSES
C
      real*8 A(500),XSQ1,XSQ2,XSQ3,XN1,XN2,XBAR1,XBAR2,TT
      INTEGER K,end1,n
        XSQ1=0.0
        XN1=0.0
        xbar1=0.0
        xbar2=0.0
       DO 20 K=1,end1
        XBAR1=XBAR1+A(K)
        XN1=XN1+1.0
 20    CONTINUE
        if (xn1.gt.0.0) then
        XBAR1=XBAR1/XN1
        else 
        xbar1=0.0
        end if
       DO 25 K=1,end1
        XSQ1=XSQ1+(A(K)-XBAR1)**2.
 25    CONTINUE
        XSQ1=XSQ1/XN1/(XN1-1.0)
        XN2=0.0
        XSQ2=0.0
       DO 30 K=end1+1,N
         XBAR2=XBAR2+A(K)
         XN2=XN2+1.0
 30    CONTINUE
         if (xn2.gt.0.0) then
         XBAR2=XBAR2/XN2
         else
         xbar2=0.0
         end if
       DO 35 K=end1+1,N
         XSQ2=XSQ2+(A(K)-XBAR2)**2.
 35    CONTINUE
         if (xn2.gt.1.0) then
         XSQ2=XSQ2/XN2/(XN2-1.0)
         XSQ3=XSQ1+XSQ2
         TT=ABS(XBAR1-XBAR2)/SQRT(XSQ3)
         else
         tt=0.0
         end if
       RETURN
       END
c
       subroutine MRPPT(num,numtot,dy,prob)
C
C         THE TEST STATISTIC ASSOCIATED WITH MULTI-RESPONSE PERMUTATION
C         PROCEDURES (MRPPT) IS CALCULATED IN THIS FORTRAN-77 PROGRAM.
C         THIS TEST IS DESIGNED TO ANALYZE MULTI-RESPONSE DATA AT
C         THE ORDINAL (OR HIGHER) LEVELS AND IS BASED ON THE AVERAGE
C         EUCLIDEAN DISTANCE BETWEEN OBJECT POSITIONS WITHIN A PRIORI
C         DISJOINT GROUPS OF OBJECTS IN A FINITE POPULATION WHERE
C         EACH POSITION IS A POINT IN AN R-DIMENSIONAL SPACE (THUS
C         THE ARGUMENTS OF EACH R-DIMENSIONAL POINT CORRESPOND TO
C         THE R MEASURED RESPONSES OF EACH OBJECT IN THE FINITE POPULATION
C         OF OBJECTS).  THIS PROGRAM ALSO CONTAINS THE TRUNCATED DISTANCE
C         OPTION (THE UPPER TRUNCATION CONSTANT IS TERMED SUP).  THE
C         PRESENT MAXIMUM VALUES OF G, N AND R IN THIS PROGRAM ARE
C         RESPECTIVELY 10, 100 AND 5.
C
C         PROGRAM MODIFIED 3/11/89
C
C       EITHER OF THE FOUR FOLLOWING WEIGHTING FORMS MAY BE CHOSEN.
C            CFORM = 1:  C(I) = n(I)/sum(n(I))
C            CFORM = 2:  C(I) = (n(I)-1)/sum(n(I)-1)
C            CFORM = 3:  C(I) = 1/sum(1)
C            CFORM = 4:  C(I) = (n(I)*(n(I)-1))/sum(n(I)*(n(I)-1))
C       CFORM IS PASSED AS THE FOURTH PARAMETER ON THE SECOND RECORD
C       IN THE DATA FILE.  WEIGHTS FOR EACH GROUP ARE CALCULATED IN
C       SUBROUTINE DIST, AND STORED IN ARRAY C(I).
C
      double precision  V,SUP,D3O,DLIT1,DLIT2,DLIT3,P1,P2,P3,P4,DBAR,
     +   C(10),XI(10),D1,VAR,T,GAM,PROB,dy(100)
      CHARACTER  HEAD*80,FIN*40
      INTEGER  CFORM,ISAM(10)
C
c      WRITE(*,*) 'ENTER FILE NAME OF INPUT DATA'
c      READ(*,1000) FIN
 1000 FORMAT(A)
c      WRITE(*,*) 'ENTER FILE NAME FOR OUTPUT RESULTS'
c      READ(*,1000) FOUT
      CALL DIST(N,ID,V,SUP,CFORM,NSAMP,ISAM,HEAD,FIN,D3O,
     +   DLIT1,DLIT2,DLIT3,P1,P2,P3,P4,DBAR,C,XI,dy,num,numtot)
      CALL CALC(N,NSAMP,ISAM,D3O,DLIT1,DLIT2,DLIT3,P1,P2,P3,P4,DBAR,
     +   C,D1,VAR,T,GAM)
c      OPEN(6,FILE=FOUT,STATUS='NEW')
c      WRITE(6,1000) HEAD
c      HEAD = 'MULTI-RESPONSE PERMUTATION PROCEDURES (MRPPT)'
c      WRITE(6,1100) HEAD
c 1100 FORMAT(/,T23,A,/)
c      WRITE(6,1200) 'INPUT HAS ', N,' OBSERVATIONS'
c 1200 FORMAT(T10,A,I3,A)
c      WRITE(6,1200) 'THERE WERE', ID, ' DIMENSION(S)'
c      WRITE(6,1300) 'IN: ', NSAMP, ' GROUPS OF SIZES ',
c     +              (ISAM(I),I = 1,NSAMP)
c 1300 FORMAT(T10,A,I5,A,10I5)
c      GO TO (1410,1420,1430,1440) CFORM
c 1410 WRITE(6,1450) 'C(I) = n(I)/sum(n(I))'
c      GO TO 1460
c 1420 WRITE(6,1450) 'C(I) = (n(I)-1)/sum(n(I)-1)'
c      GO TO 1460
c 1430 WRITE(6,1450) 'C(I) = 1/sum(1)'
c      GO TO 1460
c 1440 WRITE(6,1450) 'C(I) = (n(I)*(n(I)-1))/sum(n(I)*(n(I)-1))'
c 1450 FORMAT(T10,A)
c 1460 WRITE(6,1500)
c 1500 FORMAT(1H )
c      DO 1700 I = 1,NSAMP
c      WRITE(6,1600) 'GROUP NUMBER ', I, ' OF SIZE ', ISAM(I),
c     +                 ' HAS AN AVERAGE DISTANCE =', XI(I)
c 1600 FORMAT(T10,A,I2,A,I3,A,G16.8)
c 1700 CONTINUE
c      WRITE(6,1500)
c      WRITE(6,1800) 'THE TEST STATISTIC IS =', T
c 1800 FORMAT(T10,A,G16.8)
c      WRITE(6,1800) 'THE DISTANCE EXPONENT =', V
c      WRITE(6,1800) 'THE UPPER BOUND VALUE =', SUP
c      WRITE(6,1800) 'THE OBSERVED DELTA IS =', DBAR
c      WRITE(6,1800) 'THE EXPECTED DELTA IS =', D1
c      WRITE(6,1800) 'THE VARIANCE OF DELTA =', VAR
c      WRITE(6,1800) 'THE SKEWNESS OF DELTA =', GAM
      CALL PVALUE(T,GAM,PROB)
c      WRITE(6,1900) 'P-VALUE OF A SMALLER OR EQUAL DELTA =', PROB
c 1900 FORMAT(/,T10,A,G16.8)
c      WRITE(6,2000) '* * *   M R P P T   E N D   * * *'
c 2000 FORMAT(/,T22,A)
c      STOP 'MRPPT RESULTS IN NAMED OUTPUT FILE'
	return
      END
C
      SUBROUTINE DIST(N,ID,V,SUP,CFORM,NSAMP,ISAM,HEAD,FIN,D3O,
     +   DLIT1,DLIT2,DLIT3,P1,P2,P3,P4,DBAR,C,XI,dy,num,numtot)
C
C         SUBROUTINE DIST CALCULATES THE DISTANCES BETWEEN THE DATA
C         POINTS, PLUS THE PARAMETERS WHICH ARE DIRECTLY DEPENDENT
C         ON THE DISTANCES.
C
      double precision  V,SUP,ZERO,ONE,DLIT1,DLIT2,DLIT3,P1,P2,P3,P4,
     +   C(10),Y,XI(10),DATA(100,6),D(100,100),DSUB1(100),DSUB2(100),
     +   FK,D3O,DBAR,dy(100)
      CHARACTER  HEAD*80,FIN*40
      INTEGER  CFORM,ISAM(10)
C
      ZERO = 0.0D0
      ONE = 1.0D0
c      OPEN(8,FILE=FIN)
c      READ(8,'(A)') HEAD
c      READ(8,*) N,ID,V,SUP,CFORM,NSAMP,(ISAM(I),I = 1,10)
	n=numtot !number in sample
	id=1 !dimensions, we 1 dimensional
	v=1 !we'll deal in absolute space not squared
	sup=200 !distance upper bounds ??? can set higher
	cform=1 !several different weighting forms, we'll use #1
	nsamp=2  !the number of groups
	isam(1)=num !the number in the first group
	isam(2)=numtot-num !the number in the second group
      L = 1
      M = ISAM(L)
      DO 600 I = 1,N
c      READ(8,*,END = 700) (DATA(I,J),J = 1,ID)
 	data(i,1)=dy(i)
  400 IF(I.LE.M) GO TO 500
      L = L+1
      IF(L.LE.NSAMP) THEN
          M = M+ISAM(L)
      ELSE
          M = N
      END IF
      GO TO 400
  500 DATA(I,ID+1) = ONE*L
  600 CONTINUE
      CLOSE(8)
      GO TO 900
  700 continue!WRITE(*,800) 'E R R O R:  END-OF-FILE ENCOUNTERED.'
  800 FORMAT(/,T3,A)
      STOP
  900 DLIT1 = ZERO
      DLIT2 = ZERO
      DLIT3 = ZERO
      DBAR = ZERO
      DO 1000 I = 1,NSAMP
      XI(I) = ZERO
 1000 CONTINUE
      DO 1100 I = 1,N
      DSUB1(I) = ZERO
      DSUB2(I) = ZERO
 1100 CONTINUE
      DO 1400 I = 1,N
      DO 1300 J = 1,N
      D(I,J) = ZERO
      DO 1200 K = 1,ID
      D(I,J) = D(I,J)+(DATA(I,K)-DATA(J,K))**2
 1200 CONTINUE
      D(I,J) = D(I,J)**(V/2)
      IF(D(I,J).GT.SUP) D(I,J) = SUP
      DLIT1 = DLIT1+D(I,J)
      DLIT2 = DLIT2+D(I,J)**2
      DLIT3 = DLIT3+D(I,J)**3
      DSUB1(I) = DSUB1(I)+D(I,J)
      DSUB2(I) = DSUB2(I)+D(I,J)**2
      IS = IDNINT(DATA(I,ID+1))
      JS = IDNINT(DATA(J,ID+1))
      IF(IS.LE.NSAMP.AND.IS.EQ.JS) XI(IS) = XI(IS)+D(I,J)
 1300 CONTINUE
 1400 CONTINUE
      GO TO (1410,1430,1450,1470) CFORM
 1410 FK = ZERO
      DO 1415 I = 1,NSAMP
        FK = FK+ISAM(I)
 1415 CONTINUE
      DO 1420 I = 1,NSAMP
        if (fk.gt.0.0) then
        C(I) = ONE*ISAM(I)/FK
        else
         c(i)=0.0
        end if
 1420 CONTINUE
      GO TO 1490
 1430 FK = ZERO
      DO 1435 I = 1,NSAMP
        FK = FK+ISAM(I)
 1435 CONTINUE
      DO 1440 I = 1,NSAMP
        if (fk-nsamp.gt.0.0) then
        C(I) = ONE*(ISAM(I)-1)/(FK-NSAMP)
        else
        c(i)=0.0
        end if
 1440 CONTINUE
      GO TO 1490
 1450 DO 1460 I = 1,NSAMP
        C(I) = ONE/NSAMP
 1460 CONTINUE
      GO TO 1490
 1470 Y = ZERO
      DO 1475 I = 1,NSAMP
        Y = Y+ONE*ISAM(I)*(ISAM(I)-1)
 1475 CONTINUE
      DO 1480 I = 1,NSAMP
        if (y.gt.0.0) then
        C(I) = ONE*ISAM(I)*(ISAM(I)-1)/Y
        else
        c(i)=0.0
        end if
 1480 CONTINUE
 1490 DO 1500 I = 1,NSAMP
      XI(I) = XI(I)/ISAM(I)/(ISAM(I)-1)
      DBAR = DBAR+C(I)*XI(I)
 1500 CONTINUE
      D3O = ZERO
      DO 1700 I = 1,N-2
      DO 1710 J = I+1,N-1
      DO 1720 K = J+1,N
      D3O = D3O+D(I,J)*D(I,K)*D(J,K)
 1720 CONTINUE
 1710 CONTINUE
 1700 CONTINUE
      if (n.gt.2) then
      D3O = D3O*6/N/(N-1)/(N-2)
      else
      D30=0
      end if
      P1 = ZERO
      P2 = ZERO
      P3 = ZERO
      P4 = ZERO
      DO 1800 J = 1,N
      P1 = P1+DSUB1(J)*DSUB1(J)
      P2 = P2+DSUB1(J)*DSUB2(J)
      P4 = P4+DSUB1(J)**3
 1800 CONTINUE
      DO 1900 I = 1,N-1
      DO 1910 J = I+1,N
      P3 = P3+DSUB1(I)*DSUB1(J)*D(I,J)
 1910 CONTINUE
 1900 CONTINUE
      RETURN
      END
C
      SUBROUTINE CALC(N,NSAMP,ISAM,D3O,DLIT1,DLIT2,DLIT3,P1,P2,P3,P4,
     +   DBAR,C,D1,VAR,T,GAM)
C
C         SUBROUTINE CALC CALCULATES ALL THE REMAINING PARAMETERS
C         TOGETHER WITH THE FIRST THREE MOMENTS OF THE SAMPLING
C         DISTRIBUTION.
C
      double precision  DLIT1,DLIT2,DLIT3,P1,P2,P3,P4,VAR,T,GAM,ZERO,
     1   ONE,COF1,COF2,COF3,COF4,COF5,D1,D2,D3,D2S,D3S,D3O,D2SS,D3SS,
     2   D3OO,D3SSS,D3OOO,D3SSSS,SU2,SU3,SU4,SU5,SU6,SU7,SU8,SU9,SQV,
     3   DBAR,C(10),SIJ,SIJ2,SIJ2Q,CJ,CJ2,CJ3
      INTEGER  ISAM(10)
C
      ZERO = 0.0D0
      ONE = 1.0D0
      COF1 = ONE*N*(N-1)
      D1 = DLIT1/COF1
      D2 = DLIT2/COF1
      D3 = DLIT3/COF1
      COF2 = COF1*(N-2)
      D2S = (P1-DLIT2)/COF2
      D3S = (P2-DLIT3)/COF2
      COF3 = COF2*(N-3)
      D2SS = (DLIT1*DLIT1-D2S*COF2*4-DLIT2*2)/COF3
      D3SS = (DLIT1*DLIT2-D3S*COF2*4-DLIT3*2)/COF3
      D3OO = (P3*2-D3S*COF2*2-D3O*COF2-DLIT3)/COF3
      D3OOO = (P4-D3S*COF2*3-DLIT3)/COF3
      COF4 = COF3*(N-4)
      D3SSS = (DLIT1*D2S*COF2-D3OO*COF3*4-D3OOO*COF3*2-
     +        D3S*COF2*4-D3O*COF2*2)/COF4
      COF5 = COF4*(N-5)
      D3SSSS = (DLIT1*D2SS*COF3-D3SSS*COF4*8-D3SS*COF3*4-
     +         D3OO*COF3*8)/COF5
      SU2 = ZERO
      SU3 = ZERO
      SU4 = ZERO
      SU5 = ZERO
      SU6 = ZERO
      SU7 = ZERO
      SU8 = ZERO
      SU9 = ZERO
      DO 1100 J = 1,NSAMP
      SIJ = ONE*ISAM(J)
      SIJ2 = SIJ*(SIJ-1)
      SIJ2Q = SIJ2*SIJ2
      CJ = C(J)
      CJ2 = CJ*CJ
      CJ3 = CJ2*CJ
      SU2 = SU2+CJ2/SIJ2
      SU3 = SU3+CJ2/SIJ
      SU4 = SU4+CJ3*4/SIJ2Q
      SU5 = SU5+CJ3*(SIJ-2)*8/SIJ2Q
      SU6 = SU6+CJ3*(SIJ-2)*(SIJ-3)*8/SIJ2Q
      SU7 = SU7+CJ2*6*(ONE-CJ+CJ*(SIJ-2)*(SIJ-3)/SIJ2)/SIJ2
      SU8 = SU8+CJ2*12*((ONE-CJ)*SIJ2*(SIJ-2)+CJ*(SIJ-2)*
     +      (SIJ-3)*(SIJ-4))/SIJ2Q
      SU9 = SU9+CJ*((ONE-CJ)*(ONE-CJ*2)+CJ*(ONE-CJ)*
     +      (SIJ-2)*(SIJ-3)*3/SIJ2+CJ2*(SIJ-2)*
     +      (SIJ-3)*(SIJ-4)*(SIJ-5)/SIJ2Q)
 1100 CONTINUE
      VAR = (SU2-ONE/COF1)*(D2-D2S*2+D2SS)*2+
     +      (SU3-ONE/N)*(D2S-D2SS)*4
      SQV = DSQRT(VAR)
      GAM = ((SU4*D3+SU5*(D3S*3+D3O)+SU6*(D3OO*3+D3OOO)+
     +      SU7*D3SS+SU8*D3SSS+SU9*D3SSSS)-
     +      D1*VAR*3-D1*D1*D1)/(VAR*SQV)
      T = (DBAR-D1)/SQV
      RETURN
      END
C
      SUBROUTINE PVALUE(T,GAM,PROB)
C
C         SUBROUTINE PVALUE CALCULATES THE PROBABILITY THAT A VALUE OF
C         T IS LESS THAN OR EQUAL TO THE OBSERVED VALUE OF T.
C
      double precision  T,GAM,PROB,ZERO,ONE,PI,R,D,F,U,A,B,G1,G2,G3,
     +   E,G,H,H0,H1,H1N,H2,H3,W,X,Y,E1,E2,E3,E4,E5
C
      ZERO = 0.0D0
      ONE = 1.0D0
      PI = DATAN(ONE)*4
      IF(DABS(GAM).LT.0.01D0) GO TO 50
      R = ONE*2/DABS(GAM)
      D = R*R
      F = R*R+10
      DO 10 I = 1,9
          D = D*(R*R+I)
   10 CONTINUE
      U = (F*2-1)*DLOG(F)/2-F+DLOG(PI*2)/2-DLOG(D)+ONE/(F*12)
     +        -ONE/(F*F*F*360)
      A = R*R-1
      B = R*R*(DLOG(R)-1)-U
      G1 = 0.045D0
      G2 = 0.09D0
      G3 = 0.015D0
      H1 = ZERO
      H2 = ZERO
      W = -1.99D0/GAM
      IF(GAM.LT.ZERO) GO TO 30
      IF(T.LT.W) THEN
          PROB = ZERO
          GO TO 60
      END IF
      X = T
      Y = T+9
      DO 20 I = 1,99
          G = A*DLOG(R+X+G1*(2*I-1))-R*(X+G1*(2*I-1))+B
          H1 = H1+E(G)
          G = A*DLOG(R+X+G2*I)-R*(X+G2*I)+B
          H2 = H2+E(G)
   20 CONTINUE
      G = A*DLOG(R+X+G1*199)-R*(X+G1*199)+B
      H1N = E(G)
      G = A*DLOG(R+X)-R*X+B
      H0 = E(G)
      G = A*DLOG(R+Y)-R*Y+B
      H3 = E(G)
      PROB = ONE-G3*(H0+(H1+H1N)*4+H2*2+H3)
      GO TO 60
   30 IF(T.GT.W) THEN
          PROB = ONE
          GO TO 60
      END IF
      X = T-9
      Y = T
      DO 40 I = 1,99
          G = A*DLOG(R-X-G1*(2*I-1))+R*(X+G1*(2*I-1))+B
          H1 = H1+E(G)
          G = A*DLOG(R-X-G2*I)+R*(X+G2*I)+B
          H2 = H2+E(G)
   40 CONTINUE
      G = A*DLOG(R-X-G1*199)+R*(X+G1*199)+B
      H1N = E(G)
      G = A*DLOG(R-X)+R*X+B
      H0 = E(G)
      G = A*DLOG(R-Y)+R*Y+B
      H3 = E(G)
      PROB = G3*(H0+(H1+H1N)*4+H2*2+H3)
      GO TO 60
   50 E1 = 0.31938153D0
      E2 = -0.356563782D0
      E3 = 1.781477937D0
      E4 = -1.821255978D0
      E5 = 1.330274429D0
      H = 0.2316419D0
      W = ONE/(H*DABS(T)+1)
      G = -T*T/2
      PROB = ((((E5*W+E4)*W+E3)*W+E2)*W+E1)*W*E(G)/DSQRT(PI*2)
      IF(T.GT.ZERO) PROB = ONE-PROB
   60 RETURN
      END
C
      DOUBLE PRECISION FUNCTION E(G)
C
      double precision G
C
      E = 0.0D0
      IF(G.GE.-650) E = DEXP(G)
      RETURN
      END
	subroutine alpha(length,falpha,falpha2,talpha,talpha2)
        real*8 falpha,falpha2,talpha,talpha2
        integer length
c
c set critical values for t and f tests
c
        if ((length.ge.11).and.(length.lt.15)) then
         falpha=5.0
         falpha2=20.0
         talpha=2.2
	 talpha2=8.6
        else if ((length.ge.15).and.(length.lt.20)) then
         falpha=3.65
         falpha2=13.9
         talpha=2.13
         talpha2=4.8
        else if ((length.ge.20).and.(length.lt.25)) then
         falpha=3.26
         falpha2=9.73
         talpha=2.09
         talpha2=4.2
        else if ((length.ge.25).and.(length.lt.30)) then
         falpha=3.09
         falpha2=8.0
         talpha=2.06
         talpha2=3.85
	else if ((length.ge.30).and.(length.lt.40)) then
         falpha=2.96
         falpha2=7.3
         talpha=2.04
        else if ((length.ge.40).and.(length.lt.60)) then
         falpha=2.80
         falpha2=6.5
         talpha=2.02
         talpha2=3.55
        else if ((length.ge.60).and.(length.lt.120)) then
         falpha=2.7
         falpha2=6.0
         talpha=2.0
         talpha2=3.46
        else if (length.ge.120) then
         falpha=2.63
         falpha2=5.6
         talpha=1.98
         talpha2=3.37
        end if
c
	return
	end

c
c  ***************************************************
c

C
      SUBROUTINE CORREL(X,Y,NN,BAD,COR,VARX,VARY,N,AVGX,AVGY,SUMX,SUMY)

      REAL*8 X(500),Y(500),N,bad,cor,varx,vary,avgx,avgy,sumx,sumy
      integer nn 

c subroutine was written by TK, given to me by DE, I've evaluated it
c compared to SAS and it gives the same value. GEE that's nice to know!
c Did you expect it not to work???? DE

      N=0.
      SUMX=0.
      SUMX2=0.
      SUMY=0.
      SUMY2=0.
      SUMXY=0.
      DO 1 I=1,NN
        IF(X(I).EQ.BAD.OR.Y(I).EQ.BAD) GO TO 1
        N=N+1.0
        SUMX=SUMX+X(I)
        SUMY=SUMY+Y(I)
        SUMX2=SUMX2+(X(I)*X(I))
        SUMY2=SUMY2+(Y(I)*Y(I))
        SUMXY=SUMXY+(X(I)*Y(I))
 1    CONTINUE
c	write(*,*)'in sub: n= ',n
      IF(N.LE.2.) THEN !must have 3 data points for corr
        COR=bad
        VARX=0.
        VARY=0.
        AVGX=0.
        AVGY=0.
        SUMY=0.
        SUMX=0.
        RETURN
      END IF
      AVGX=SUMX/N
      AVGY=SUMY/N
      VARX=(SUMX2-(N*((SUMX/N)**2.)))/(N-1.)
      VARY=(SUMY2-(N*((SUMY/N)**2.)))/(N-1.)
      COR=(N*SUMXY - SUMX*SUMY)/SQRT((N*SUMX2 - SUMX*SUMX)*
     1  (N*SUMY2 - SUMY*SUMY))

 98	format(3f8.4)

      RETURN
 99   END
c***********************************************************************sta


c     call good_nuf(miss,a,b,igoodnuf,nummy)
      subroutine good_nuf(zindef,cand,ref,igoodnuf,nummy)

c  this subroutine takes the annual reference series output,
c     handed down to it
c  checks to see if:
c     a.  there are 10 or more years with good data & reference series
c     b.  the correlation of the two is .ge. 0.8  which is our cut off
c  (see comments in ep.f) for when we can evaluate these time series

c     returns igoodnuf:  0 equals OK, 1 equals not good enough
c     also shortens the reference series array to see if early years
c     were the ones causing a problem

      integer igoodnuf
      real*8 cand(500),ref(500)
      real*8 r(2,500),x(500),y(500)
      real*8 zindef !missing value


c     iy1&iy2 are start and end years

      zin=zindef+1


 50   continue

      iy1=9999
         do j=1,nummy
            r(1,j)=cand(j)
            r(2,j)=ref(j)
            if(cand(j).gt.zin.and.ref(j).gt.zin)then
               iy1=min(iy1,j)
               iy2=j
            endif
        enddo

   
 100  continue


c     enough data? 


c     20 years of data required:
      num=0
      do i=iy1,iy2
         if(r(1,i).gt.zin.and.r(2,i).gt.zin)num=num+1
      enddo
      do i=iy1,iy2
         if(r(1,i).gt.zin.and.r(2,i).gt.zin)then
            ky1=i !first year of data for both
            go to 287
         endif
      enddo
 287  continue
      do i=iy2,iy1,-1
         if(r(1,i).gt.zin.and.r(2,i).gt.zin)then
            ky2=i !last year of data
            go to 288
         endif
      enddo
 288  continue
      numk=1+ky2-ky1 !first to last year must be 20 years 
      if(numk.lt.20.or.num.lt.15)then !because one month might be missing
c in the middle somewhere which would be ok, 
c we limit it to 20 yrs with some missing and
c no more than 5 years of that missing final annual values
         igoodnuf=1 !not good enough
         go to 500
      endif
c     needs at least 10 years of data and reference series to make
c     some adjustments
      num=0
      do i=iy1,iy2
         if(r(1,i).gt.zin.and.r(2,i).gt.zin)num=num+1
      enddo
      if(num.lt.10)then
         igoodnuf=1 !not good enough
         go to 500
      endif
c ok, enough data and reference series.  But are the reference series
c good enough.  The cut off value is 0.80 (rounded up from 0.79, so we
c will use .lt.0.8)
      do i=1,500
         x(i)=zindef
         y(i)=zindef
      enddo
      num=0
      do i=iy1,iy2-1
         if(r(1,i).gt.zin.and.r(2,i).gt.zin)then
            if(r(1,i+1).gt.zin.and.r(2,i+1).gt.zin)then
               num=num+1
               x(num)=r(1,i+1)-r(1,i)  !dT/dt series correlations
               y(num)=r(2,i+1)-r(2,i)
            endif
         endif
      enddo

       call CORREL8(X,Y,num,zindef,COR,VARX,VARY,aN,AVGX,AVGY,SUMX,SUMY)
	
      if(cor.lt.0.80)then
c   presumably the early years would have worse reference series
c   because of less station density.  So, we'll remove the earliest
c   years one at a time until we get one that works or it is too 
c   short to 
         ref(iy1)=zindef
         go to 50
      endif

      if(cor.ge.0.80)then
         igoodnuf=0
         go to 500
      endif

 500  continue
      return
      end

************************************************




      SUBROUTINE CORREL8
     +(X,Y,NN,zindef,COR,VARX,VARY,aN,AVGX,AVGY,SUMX,SUMY)

      integer NN
      REAL*8 X(500),Y(500),N 
c note, N is real not integer
c	NN is the size of the array
	N=aN ! changed in call statement to ensure it not being integer
	bad=zindef !missing value passed down, assumed to be very negative
c subroutine was written by TK, given to me by DE, I've evaluated it
c compared to SAS and it gives the same value



      N=0.
      SUMX=0.
      SUMX2=0.
      SUMY=0.
      SUMY2=0.
      SUMXY=0.
      DO 1 I=1,NN
        IF(X(I).lt.BAD+1.OR.Y(I).lt.BAD+1) GO TO 1
        N=N+1
        SUMX=SUMX+X(I)
        SUMY=SUMY+Y(I)
        SUMX2=SUMX2+(X(I)*X(I))
        SUMY2=SUMY2+(Y(I)*Y(I))
        SUMXY=SUMXY+(X(I)*Y(I))
 1    CONTINUE
      IF(N.LE.2.) THEN !must have 3 data points for corr
        COR=bad
        VARX=0.
        VARY=0.
        AVGX=0.
        AVGY=0.
        SUMY=0.
        SUMX=0.
        RETURN
      END IF
      AVGX=SUMX/N
      AVGY=SUMY/N
      VARX=(SUMX2-(N*((SUMX/N)**2.)))/(N-1.)
      VARY=(SUMY2-(N*((SUMY/N)**2.)))/(N-1.)
      COR=(N*SUMXY - SUMX*SUMY)/SQRT((N*SUMX2 - SUMX*SUMX)*
     1  (N*SUMY2 - SUMY*SUMY))
      aN=N
      RETURN
 99   END









