! scalemod.f90
! module procedure written by Tim Mitchell in Jan 2000
! last modification on 10.02.00
! includes all the routines necessary to run a scaling operation
!	PredictRun, PredictError, PredictErrorStat

module ScaleMod

implicit none

contains

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

subroutine PredictRun (JobMemN,JobYearN,ActualExe,PredictWye,Aye)

real, dimension (:,:), pointer 	:: ActualExe, PredictWye
real, intent (in) 		:: Aye
integer, intent (in) 		:: JobMemN, JobYearN

real, parameter :: MissVal = -999.0

integer :: XYear, XMem

do XMem = 1, JobMemN
  do XYear = 1, JobYearN
    if (ActualExe(XMem,XYear).NE.MissVal) then
      PredictWye(XMem,XYear) = Aye * ActualExe(XMem,XYear)
    end if
  end do
end do

end subroutine PredictRun

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

subroutine PredictError (MissAccept,RefMethod,JobMemN,JobYearN,SliceLen,ActualWye,PredictWye,&
			 ErrorMSE,ReferenceMSE,FractionMSE)

real, dimension (:,:), pointer 	:: ActualWye, PredictWye
real, dimension (:)  , pointer  :: ErrorMSE, ReferenceMSE, FractionMSE

real, intent (in)		:: MissAccept

integer, intent (in) 		:: JobMemN, JobYearN, SliceLen
integer, intent (in) 		:: RefMethod	! 1=internal-var 2=controlzero 3=Brier

real, dimension (:,:), allocatable 	:: ActualMean, PredictMean

integer, dimension (:,:), allocatable	:: YearCheck
integer, dimension (:), allocatable	:: MemCheck

real, parameter :: MissVal = -999.0

real :: ActualTot, PredictTot, DiffSqTot, ValidEn

integer :: XMem, XMem0, XMem1, XYear, XPredict, XActual
integer :: SliceMidYear, SliceStartYear, SliceEndYear
integer :: Threshold
integer :: TotalMiss, TotalValid
integer :: AllocStat

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

ErrorMSE     = MissVal
ReferenceMSE = MissVal
FractionMSE  = MissVal

allocate (ActualMean (JobMemN,JobYearN), &
	  PredictMean(JobMemN,JobYearN), &
	  YearCheck  (JobMemN,JobYearN), &
	  MemCheck   (JobYearN),         stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

ActualMean  = MissVal
PredictMean = MissVal
YearCheck   = 0
MemCheck    = 0

Threshold = nint (MissAccept * real(SliceLen)/ 100.0)

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

do XMem = 1, JobMemN
  do SliceStartYear = 1, (JobYearN-SliceLen+1)
    SliceMidYear = SliceStartYear + floor(real(SliceLen)/2.0)
    TotalMiss    = 0
    ValidEn      = 0.0
    ActualTot    = 0.0
    PredictTot   = 0.0
    
    do XYear = SliceStartYear, (SliceStartYear+SliceLen-1)
      if (ActualWye(XMem,XYear).EQ.MissVal.OR.PredictWye(XMem,XYear).EQ.MissVal) then
      	TotalMiss  = TotalMiss  + 1
      else
        ActualTot  = ActualTot  + ActualWye(XMem,XYear)
        PredictTot = PredictTot + PredictWye(XMem,XYear)
        ValidEn    = ValidEn    + 1.0
      end if
    end do
    
    if (TotalMiss.LE.Threshold.AND.ValidEn.GE.1) then	! checks that no. of valid values in slice >= 90%
    	YearCheck   (XMem,SliceMidYear) = 1
    	ActualMean  (XMem,SliceMidYear) = ActualTot  / ValidEn
    	PredictMean (XMem,SliceMidYear) = PredictTot / ValidEn
    end if
  end do
end do

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

do SliceStartYear = 1, (JobYearN-SliceLen+1)
    SliceMidYear = SliceStartYear + floor(real(SliceLen)/2.0)
    TotalValid   = 0
    
    do XMem = 1, JobMemN
      if (YearCheck(XMem,SliceMidYear).EQ.1) TotalValid = TotalValid + 1
    end do
    
    if (TotalValid.GT.1) MemCheck (SliceMidYear) = 1
end do

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

do SliceStartYear = 1, (JobYearN-SliceLen+1)
  SliceMidYear = SliceStartYear + floor(real(SliceLen)/2.0)
  
  if (MemCheck(SliceMidYear).EQ.1) then		! if sufficient members for a calc
    DiffSqTot    = 0.0
    ValidEn      = 0.0
    do XPredict = 1, JobMemN			! error
     do XActual = 1, JobMemN
      if (YearCheck(XMem,SliceMidYear).EQ.1) then
       if (PredictMean(XPredict,SliceMidYear).NE.MissVal.AND.ActualMean(XActual,SliceMidYear).NE.MissVal) then
        DiffSqTot = DiffSqTot + (PredictMean(XPredict,SliceMidYear) - ActualMean(XActual,SliceMidYear)) ** 2
        ValidEn   = ValidEn   + 1.0
       end if
      end if
     end do
    end do
    
    if (ValidEn.GT.1) ErrorMSE (SliceMidYear) = DiffSqTot / ValidEn
    
    if (RefMethod.EQ.1) then			! internal var = reference
    
     DiffSqTot    = 0.0
     ValidEn      = 0.0
     do XPredict = 1, JobMemN
      do XActual = 1, JobMemN
       if (YearCheck(XMem,SliceMidYear).EQ.1) then
        if (ActualMean(XActual,SliceMidYear).NE.MissVal) then
         DiffSqTot = DiffSqTot + (ActualMean(XPredict,SliceMidYear) - ActualMean(XActual,SliceMidYear)) ** 2
         ValidEn   = ValidEn   + 1.0
        end if
       end if
      end do
     end do
    
     if (ValidEn.GT.JobMemN) ReferenceMSE (SliceMidYear) = DiffSqTot / ValidEn
    
    else if (RefMethod.EQ.2.OR.RefMethod.EQ.3) then		! control of pred=0 = reference
    
     DiffSqTot    = 0.0
     ValidEn      = 0.0
     do XActual = 1, JobMemN
       if (YearCheck(XMem,SliceMidYear).EQ.1) then
        if (ActualMean(XActual,SliceMidYear).NE.MissVal) then
         DiffSqTot = DiffSqTot + (ActualMean(XActual,SliceMidYear)) ** 2
         ValidEn   = ValidEn   + 1.0
        end if
       end if
     end do
    
     if (ValidEn.GT.1) ReferenceMSE (SliceMidYear) = DiffSqTot / ValidEn
    
    end if 
    						! calc of fraction
    
    if (ReferenceMSE(SliceMidYear).NE.MissVal.AND.ErrorMSE(SliceMidYear).NE.MissVal) then
     if (ReferenceMSE(SliceMidYear).NE.0) then
      if (RefMethod.EQ.1.OR.RefMethod.EQ.2) then
        FractionMSE (SliceMidYear) = ErrorMSE (SliceMidYear) / ReferenceMSE (SliceMidYear)
      else if (RefMethod.EQ.3) then
        FractionMSE (SliceMidYear) = 1 - (ErrorMSE (SliceMidYear) / ReferenceMSE (SliceMidYear))    
      end if
     end if
    end if
    
  end if
  
end do

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

deallocate (ActualMean,PredictMean,YearCheck,MemCheck)

end subroutine PredictError

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

subroutine SetUpErrorStat (MissAccept,JobADYear,JobYearN,SliceLen,QMethod,FoundYear0,FoundYear1,&
			     StatThreshold,MissThreshold)

integer, dimension (:), pointer :: JobADYear

real, intent (in)     :: MissAccept
real, intent (out)    :: StatThreshold
integer, intent (in)  :: JobYearN, SliceLen
integer, intent (out) :: FoundYear0,FoundYear1,MissThreshold
integer, intent (out) :: QMethod	! 1=%>threshold 2=mean-over-period

real, parameter :: MissVal = -999.0

integer	:: LowerTrim, UpperTrim		! allowances in MissThreshold due to trimming from SliceLen>1 
integer :: ReadStatus
integer :: GetYear0,GetYear1
integer :: XYear

!***************************************
! make selections

!print*, "  > Select the prediction error statistic to calculate: "
!print*, "  >   1. % of chosen period that Error/Internal < chosen threshold"
!print*, "  >   2. mean of Error/Internal over chosen period"
!do
!	read (*,*,iostat=ReadStatus), QMethod
!	if (ReadStatus.LE.0.AND.QMethod.GE.1.AND.QMethod.LE.2) exit
!end do

QMethod = 2

print*, "  > Select first and last years AD of period: "
do
	read (*,*,iostat=ReadStatus), GetYear0, GetYear1
	if (ReadStatus.LE.0.AND.GetYear1.GE.GetYear0) exit
end do

if (QMethod.EQ.1) then
  print*, "  > Select the Error/Internal threshold : "
  do
	read (*,*,iostat=ReadStatus), StatThreshold
	if (ReadStatus.LE.0) exit
  end do
else
  StatThreshold = MissVal
end if

!***************************************
! locate years AD

XYear = 0
FoundYear0 = 0
do
  XYear = XYear + 1
  
  if (JobADYear(XYear).EQ.GetYear0) FoundYear0 = XYear
  
  if (FoundYear0.GT.0)   exit
  if (XYear.EQ.JobYearN) exit
end do
if (FoundYear0.EQ.0) print*, "  > ##### ERROR: Initial year selected not found #####"

FoundYear1 = FoundYear0 + GetYear1 - GetYear0

UpperTrim = max ((FoundYear1-(JobYearN-floor(real(SliceLen/2.0)))),0)
LowerTrim = max ((floor(real(SliceLen/2.0))-FoundYear0),0)

MissThreshold = nint (MissAccept * (real(FoundYear1 - FoundYear0 + 1) / 100.0)+UpperTrim+LowerTrim)

end subroutine SetUpErrorStat

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

subroutine PredictErrorStat (AllFractionMSE,JobADYear,RegN,JobYearN,QMethod,FoundYear0,FoundYear1,&
			     StatThreshold,MissThreshold,AllErrorStat)

real, dimension (:,:), pointer 	:: AllFractionMSE
real, dimension (:), pointer	:: AllErrorStat
integer, dimension (:), pointer :: JobADYear

real, intent (in)		:: StatThreshold
integer, intent (in) 		:: RegN, JobYearN, QMethod, FoundYear0, FoundYear1, MissThreshold

real, parameter :: MissVal = -999.0

real :: PerEn, PerTot

integer :: XYear, XReg
integer :: ReadStatus
integer :: PerMiss

!***************************************
! method 1

if (QMethod.EQ.1) then
  do XReg = 1, RegN
    
    PerTot  = 0.0
    PerEn   = 0.0
    PerMiss = 0
  
    do XYear = FoundYear0, FoundYear1
      if (AllFractionMSE(XReg,XYear).NE.MissVal) then
        PerEn   = PerEn  + 1.0
        
        if (AllFractionMSE(XReg,XYear).LT.StatThreshold) &
          	PerTot  = PerTot + 1.0
      else
        PerMiss = PerMiss + 1
      end if
    end do
    
    if (PerMiss.LE.MissThreshold) then
      AllErrorStat (XReg) = 100.0 * (PerTot / PerEn)
    else
      AllErrorStat (XReg) = MissVal
    end if

  end do
    
end if

!***************************************
! method 2

if (QMethod.EQ.2) then
  do XReg = 1, RegN
    
    PerTot  = 0.0
    PerEn   = 0.0
    PerMiss = 0
  
    do XYear = FoundYear0, FoundYear1
      if (AllFractionMSE(XReg,XYear).NE.MissVal) then
        PerTot  = PerTot + AllFractionMSE(XReg,XYear)
        PerEn   = PerEn  + 1.0
      else
        PerMiss = PerMiss + 1
      end if
    end do
    
    if (PerMiss.LE.MissThreshold) then
      AllErrorStat (XReg) = PerTot / PerEn
    else
      AllErrorStat (XReg) = MissVal
    end if

  end do
    
end if

end subroutine PredictErrorStat

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

end module ScaleMod