        
*	************************************************
*	
*	       STATISTICAL FORTRAN SUBROUTINES 
*	
*	          Byron Gleason, 4/2000
*
*	************************************************

*
*	**************
*
*	INDEX OF STATISTICAL ROUTINES WITH BRIEF DESCRIPTION
*	(For more detailed information skip down to the actual routine)
*
*
*	1. LSQREG = Calculates a linear least squares regression and returns
*	 	    a slope and intercept.
*
*	2. TREND - Calculates a least squares regression and then finds
*	  	   the predicted y-value for a given x-value.
*
*	3. AVERAGE - Calculates the average value of an array of reals
*
*
*	**************


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

	subroutine LSQREG(n,miss,xdata,ydata,slope,yint)

 	integer n 
	real sumx,sumx2,sumy,sumxy,xmean,ymean,slope,yint,
     +  xdata(1:n),ydata(1:n),miss

	count = 0
	sumx = 0
	sumx2 = 0 
	sumy = 0
	sumxy = 0
	
	do i=1,n
	  if ((xdata(i).eq.miss).or.(ydata(i).eq.miss)) then
	  else
	    count = count + 1.0
	    sumx = sumx + xdata(i)
	    sumx2 = sumx2 + xdata(i) ** 2
	    sumy = sumy + ydata(i)
	    sumxy = sumxy + xdata(i) * ydata(i) 
	  endif
	enddo	  

c	calculate slope and y-intercept

	if (count.eq.0) then
	  slope=miss
	  yint=miss
	else
	  xmean = sumx/count
	  ymean = sumy/count
	  denom = (sumx2-sumx*xmean)
	  if (denom.ne.0) then
	    slope = (sumxy-sumx*ymean) / (sumx2-sumx*xmean)
	    yint = ymean - slope * xmean
	  else
	    slope = miss
	    yint = miss
	  endif
	endif
	return
	end

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

	subroutine TREND(n,rMiss,xdata,ydata,xvalue,yvalue)

	real xdata(1:n),ydata(1:n),xvalue,yvalue,slope,yint,
     +  rMiss


	call lsqreg(n,rMiss,xdata,ydata,slope,yint)

c	calculate yvalue using equation of the form: Y = MX + B
c	where M = slope, B = y-intercept (yint), xvalue = X

	if (slope.eq.rMiss) then
	  yvalue=rMiss
	else
	  yvalue = (slope*xvalue) + yint
	endif

	return
	end

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

	subroutine AVERAGE(n,rMiss,rData,rAve)

	real rData(1:n)

	rCount=0
	rSum = 0.0
	
	do i=1,n
	  if (rData(i).ne.rMiss) then
	    rCount=rCount+1.
	    rSum = rSum+rData(i)
	  endif
	enddo

	if (rCount.gt.0.) then
	  rAve=(rSum/rCount)
	else
	  rAve=rMiss
	endif

	return
	end

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