! compare.f90
! written by Tim Mitchell on 20.03.03
! module contains routines to clean-up .dtb files through comparisons

module Compare

use FileNames
use PlainPerm
use Regress

implicit none

contains

!*******************************************************************************
! produces monthly vectors (a,b,weight) for StnX and StnY, under y=a+bx
! beg,end specify the limits of the year subscripts to use
! if MonthR is set, Pearson's r is returned as well as weight
! if MeanR is set, an annual 'r' estimate is returned; and if ThreshR is set, 
!	every month must be present and with 'r' > ThreshR, or else = MissVal

subroutine PairRelation (Differ,StnX,StnY,MonthA,MonthB,MonthR, &
			 MonthSX,Beg,End,MeanR,ThreshR)

integer, dimension (:,:), pointer 		:: StnX,StnY
real, dimension (:), pointer 			:: Exe,Wye
real, dimension (12), intent(out) 		:: MonthA,MonthB,MonthR
real, dimension (12), intent(out), optional 	:: MonthSX	! pop stdev of X
logical, intent(in)				:: Differ
real, intent(out), optional			:: MeanR
real, intent(in), optional			:: ThreshR
integer, intent(in), optional			:: Beg,End

integer, parameter :: DataMissVal=-9999
real, parameter :: MissVal=-999.0
integer, parameter :: NMonth=12

real    :: Are,Enn, OpTot,OpTotSq,OpEn,OpVal,CalcThreshR
integer :: XMonth,XYearIn,NYearIn,XYearSub,NYearSub,Year0,Year1
integer :: AllocStat,QValid,Elastic

Elastic=10 ; if (Differ) Elastic=5
NYearIn=size(StnX,1) ; QValid=1
if (NYearIn.NE.size(StnY,1)) print*, "  > ##### ERROR: PairRelation size"

Year0=1 ; Year1=NYearIn
if (present(Beg).AND.present(End)) then
  Year0=Beg ; Year1=End
  if (Year0.LT.1.OR.Year1.GT.NYearIn.OR.Year0.GT.Year1) &
  		print*, "  > ##### ERROR: PairRelation range", Year0,Year1
end if
NYearSub=Year1-Year0+1

allocate (Exe(NYearSub),Wye(NYearSub),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PairRelation: Allocation failure #####"
MonthR = MissVal

CalcThreshR=MissVal-1

QValid=1 ; XMonth=0
do
  XMonth=XMonth+1
  
  Exe=MissVal ; Wye=MissVal
  do XYearIn = Year0,Year1
    XYearSub = XYearIn-Year0+1
    if (StnX(XYearIn,XMonth).NE.DataMissVal) Exe(XYearSub)=StnX(XYearIn,XMonth) 
    if (StnY(XYearIn,XMonth).NE.DataMissVal) Wye(XYearSub)=StnY(XYearIn,XMonth)
  end do

  call RegressVectors (Exe,Wye,MonthA(XMonth),MonthB(XMonth),MonthR(XMonth),Enn,Elastic)
  
  if (MonthR(XMonth).LE.CalcThreshR) QValid=0
  if (XMonth.EQ.12.OR.QValid.EQ.0) exit
end do

if (present(MeanR)) then
  if (QValid.EQ.1) then
    MeanR=0
    do XMonth = 1, NMonth
      MeanR=MeanR+MonthR(XMonth)
    end do
    MeanR=MeanR/12.0
  else
    MeanR = MissVal
  end if
end if

deallocate (Exe,Wye,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PairRelation: Deallocation failure #####"

end subroutine PairRelation

!*******************************************************************************

end module Compare