subroutine ftest(data1,n1,data2,n2,f,prob)
      implicit none
      integer n1,n2
      double precision f,prob,data1(n1),data2(n2)
!given the arrays data1 and data2, this routine returns the value of f
!and the significance as prob. Small values of prob indicate that the
!two arrays have significantly different variances
      double precision ave1,ave2,df1,df2,var1,var2,betai,
     &xmin,xmax
      call stdmaxmeanmin(n1,data1,var1,ave1,xmin,xmax)
      var1=var1*var1
      call stdmaxmeanmin(n2,data2,var2,ave2,xmin,xmax)
      var2=var2*var2
      if(var1.gt.var2)then
        f=var1/var2
        df1=dble(n1-1)
        df2=dble(n2-1)
      else
        f=var2/var1
        df1=dble(n2-1)
        df2=dble(n1-1)
      endif
      prob=2.0d0*betai(0.5d0*df2,0.5d0*df1,df2/(df2+df1*f))
      if(prob.gt.1.0d0)prob=2.0d0-prob
      return
      end

c####################################################################
      subroutine ftestrsq_rms(ymeas,ypred,n0,nparams,rsq,rms,agrind,
     &rmse_norm,rmse_perc,aic,aicc,f,prob)
      implicit double precision (a-h,l,o-z)
      dimension ymeas(n0),ypred(n0),y1(n0),y2(n0)
      fn9999=-9999.0d0
      tiny=1.0d-7
      n=0
      do i=1,n0
        if(dabs(ymeas(i)-fn9999).gt.tiny.and.
     &dabs(ypred(i)-fn9999).gt.tiny)then
          n=n+1
          y1(n)=ymeas(i)
          y2(n)=ypred(i)
        endif
      enddo
      ymin=y1(1)
      ymax=y1(1)
      do i=2,n
        if(y1(i).lt.ymin)ymin=y1(i)
        if(y1(i).gt.ymax)ymax=y1(i)
      enddo
      sum=0.0d0
      rmse_perc=0.0d0
      do 10 i=1,n
        sum=sum+(y1(i)-y2(i))*(y1(i)-y2(i))
        rmse_perc=rmse_perc+(y1(i)-y2(i))*(y1(i)-y2(i))/(y2(i)*y2(i))
10    continue   
      rms=dsqrt(sum/dble(n))
      if(nparams.gt.0)then
        aic=dble(n)*dlog(rms*rms)+2.0d0*dble(nparams)
        aicc=aic+2.0d0*dble(nparams*(nparams+1))/dble(n-nparams-1)
      else
        aic=-9999.0d0
        aicc=-9999.0d0
      endif
      rmse_norm=rms/(ymax-ymin)
      rmse_perc=100.0d0*dsqrt(rmse_perc/dble(n))
      ymean1=0.0d0
      ymean2=0.0d0
      do 20 i=1,n
        ymean1=ymean1+y1(i)
        ymean2=ymean2+y2(i)
20    continue
      ymean1=ymean1/dble(n)
      ymean2=ymean2/dble(n)
      sum1=0.0d0
      sum2=0.0d0
      sum3=0.0d0
      sum4=0.0d0
      sum5=0.0d0
      do 30 i=1,n
        sum1=(y1(i)-ymean1)*(y2(i)-ymean2)+sum1
        sum2=(y1(i)-ymean1)*(y1(i)-ymean1)+sum2
        sum3=(y2(i)-ymean2)*(y2(i)-ymean2)+sum3
        sum4=(y1(i)-y2(i))*(y1(i)-y2(i))+sum4
        sum5=(dabs(y2(i)-ymean1)+dabs(y1(i)-ymean1))*
     &(dabs(y2(i)-ymean1)+dabs(y1(i)-ymean1))+sum5
        sum6=sum6+(y2(i)-ymean1)*(y2(i)-ymean1)
30    continue
      if((sum2*sum3).eq.0.0d0)then
        rsq=-9999.0d0
      else
        rsq=sum1/dsqrt(sum2*sum3)
        rsq=rsq*rsq
      endif
      if(sum5.eq.0.0d0)then
        agrind=-9999.0d0
      else
        agrind=1.0d0-sum4/sum5
      endif
      if(rsq.gt.0.0d0)then
        df1=dble(nparams-1)
        df2=dble(n-nparams)
        f=(sum6/df1)/(sum4/df2)
        prob=betai(0.5d0*df2,0.5d0*df1,df2/(df2+df1*f))
      else
        f=-9999.0d0
        prob=-9999.0d0
      endif
      return
      end
!####################################################################