! scaleops.f90
! f90 main program written on 10.01.00 by Tim Mitchell
! last modification on 09.05.00
! f90 -o scaleops initialmod.f90 loadmod.f90 scalemod.f90 savemod.f90 transformmod.f90 scaleops.f90

program ScaleOps

use InitialMod
use LoadMod
use ScaleMod
use SaveMod
use TransformMod

implicit none

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

real, pointer, dimension (:)		:: WorkRegAye, WorkAllErrorStat
real, pointer, dimension (:)		:: TempUnSmoo, TempLowSmoo, TempHighSmoo
real, pointer, dimension (:)		:: WorkErrorMSE, WorkInternalMSE, WorkFractionMSE
real, pointer, dimension (:)		:: WorkDiffErrorMSE, WorkDiffInternalMSE, WorkDiffFractionMSE
real, pointer, dimension (:)		:: WorkSpaceA, WorkSpaceB, WorkSpaceR
real, pointer, dimension (:)		:: WorkInExe, WorkInWye, WorkOutReg, RegToBeSaved
real, pointer, dimension (:,:)		:: TimSeries, WorkAllFractionMSE, LinToBeSaved 
real, pointer, dimension (:,:)		:: WorkTimeABR, WorkSpaceABR 
real, pointer, dimension (:,:)		:: WorkActualExe, WorkActualWye, WorkPredictWye, WorkPredictor
real, pointer, dimension (:,:)		:: WorkOutExe, WorkOutWye, WorkOutTim
real, pointer, dimension (:,:,:)	:: WorkPredictand

integer, pointer, dimension (:)		:: WorkADYear, WorkMapRawReg, WorkRegSizes
integer, pointer, dimension (:,:)	:: WorkMapIDLRaw, WorkMapIDLReg

character (len=20), pointer, dimension (:) 	:: WorkRegNames 
character (len=80), pointer, dimension (:) 	:: TimRegNames, LinToBeSavedNames

real, allocatable, dimension (:,:,:)	:: WorkAllPredict, WorkAllActual

integer, allocatable, dimension (:) 	:: SaveRegKey, MemKey

character (len=10), allocatable, dimension (:)	:: MemberName

real, parameter :: MissVal = -999.0

real :: WorkAye, WorkBee, WorkPea, ErrStatStatThreshold
real :: LSRAye, LSRBee, LSRPea
real :: RegTotal, RegEn, YearTotal, YearEn, Threshold, ActualTotal, PredictTotal, MemEn
real :: WorkMissAccept		! the % acceptable missing values ; default = 10
real :: ParticAye
real :: OpTotal, OpEn
real :: MemThresh

integer :: WorkGrid,WorkLongN,WorkLatN,WorkDataN
integer :: WorkYearN,WorkDecN
integer :: WorkRegN
integer :: SaveLinN
integer :: XYear, XMem, XReg, XLine, XDec, XPredictand, XPredictor, XGloFile
integer :: ReadStatus, AllocStat
integer :: WorkMemN
integer :: TimRegN
integer :: PredictorReg
integer :: MainMenuChoice
integer :: GaussChoice
integer :: SaveRegSelect
integer :: TotRegSelect
integer :: PredictActive	! 1=Predictor and Predictand are activated 
integer :: AyeActive		! 1=WorkRegAye is activated
integer :: ErrStatMethod,ErrStatYear0,ErrStatYear1, ErrStatMissThreshold
integer :: WorkSliceLen
integer :: SaveErrStat, SaveCorrStat, SaveDiffStat
integer :: SaveTimeABR, SaveSpaceA, SaveSpaceB, SaveSpaceR
integer :: SelectMem, SelectMemN, XSelectMem
integer :: SelectADYear, SelectXYear, SelectDiffLin, SelectNatLin
integer :: QEnsSlice, QDiffScat, QGloRMSE, QDiffTim, QIntraVar, QSampPop
integer :: TotalMiss, TotalMissAye, TotalMissBee, TotalMissPea

character (len=10) :: WorkGridTitle
character (len=40) :: WorkRegTitle
character (len=80) :: WorkGridFilePath, WorkGloTitle, SaveLinTitle

!*******************************************************************************
! main routine

open (99,file="/cru/u2/f709762/data/scratch/log-scale.txt",status="replace",action="write")

call Preliminaries

do
  call MainMenu
  
  if (MainMenuChoice.EQ. 1) then
    call EnsembleSize
    call GetTimFiles
  end if
  
  if (MainMenuChoice.EQ. 2) call CalcAye
  if (MainMenuChoice.EQ. 3) call SaveAyeGlo
  if (MainMenuChoice.EQ. 4) call LoadAyeGlo  
  if (MainMenuChoice.EQ. 5) call ScalingPrediction
  if (MainMenuChoice.EQ. 6) call CalcNatVar
  if (MainMenuChoice.EQ. 7) call SetMissAccept (WorkMissAccept)
  
  if (MainMenuChoice.EQ.99) exit
end do

call WindDown

close (99)

contains

!*******************************************************************************
! main menu

subroutine MainMenu

print*, "  > Make your selection: "
print*, "  >   1. Load ensemble of model predictors and predictands."
print*, "  >   2. Use current ensemble to construct scaling equation by LSR."
print*, "  >   3. Save scaling equation to .glo file."
print*, "  >   4. Load scaling equation from .glo file."
print*, "  >   5. Apply current scaling equation to current ensemble."
print*, "  >   6. Examine nat var of current ensemble."
print*, "  >   7. Alter acceptable % of MissVal, currently: ", WorkMissAccept
print*, "  >  99. Exit."

do
	read (*,*,iostat=ReadStatus), MainMenuChoice
	if (ReadStatus.LE.0 .AND. MainMenuChoice.GE.1) exit
end do

end subroutine MainMenu

!*******************************************************************************
! preliminaries

subroutine Preliminaries

WorkMissAccept = 10.0		! default value for acceptable missing values

print*
print*, "  > ***** ScaleOps *****"
print*, "  > Performs scaling operations"
print*

print*, "  > Make initial selections that will govern all operations."

call GridSelect   (WorkGrid,WorkGridTitle,WorkLongN,WorkLatN,WorkDataN,WorkGridFilePath)
call PeriodSelect (WorkYearN,WorkDecN,WorkADYear)
call RegSelect    (WorkGrid,WorkLongN,WorkLatN,WorkDataN,WorkMapIDLReg,WorkRegSizes,WorkRegNames,&
		   WorkRegTitle,WorkRegN)

print*

end subroutine Preliminaries

!*******************************************************************************
! specify ensemble size

subroutine EnsembleSize

if (PredictActive.EQ.1) then
  deallocate (WorkPredictor,WorkPredictand,stat=AllocStat)
  if (AllocStat.EQ.0) PredictActive = 0
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure #####"
end if

print*, "  > Enter the number of runs to load: "
do
	read (*,*,iostat=ReadStatus), WorkMemN
	if (ReadStatus.LE.0 .AND. WorkMemN.GE.1) exit
end do

allocate ( WorkPredictor  (WorkMemN,WorkYearN), &
	   WorkPredictand (WorkRegN,WorkMemN,WorkYearN), stat=AllocStat)
if (AllocStat.EQ.0) PredictActive = 1
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

WorkPredictor  = MissVal
WorkPredictand = MissVal

allocate ( MemberName  (WorkMemN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

end subroutine EnsembleSize

!*******************************************************************************
! get .tim files

subroutine GetTimFiles

print*, "  > Load time series from .tim files."
print*, "  > Predictor is independent of choices made thus far, except years."
print*, "  > Predictand must be of correct region set and years."

do XMem = 1, WorkMemN
  print*, "  > *** MEMBER ", XMem
  print*, "  > Name this member: "
  
  do
	read (*,*,iostat=ReadStatus), MemberName(XMem)
	if (ReadStatus.LE.0.AND.MemberName(XMem).NE."") exit
  end do
  
  write (99,"(a8,a20)"), "member: ", MemberName(XMem)

  print*, "  > PREDICTOR:"
  call LoadTim (WorkYearN,WorkADYear,TimRegN,TimRegNames,TimSeries)
    
  if (TimRegN.GT.1) then
    	print*, "  > Select a region (0=display list): "
    	do
		read (*,*,iostat=ReadStatus), PredictorReg
		if (PredictorReg.EQ.0) then
		  do XReg = 1, TimRegN
		    print "(I6,A1,A40)", XReg, " ", trim(adjustl(TimRegNames(XReg)))
		  end do
		end if
		if (ReadStatus.LE.0.AND.PredictorReg.GE.1.AND.PredictorReg.LE.TimRegN) exit
    	end do
  else
    	PredictorReg = 1
  end if
  
  TotalMiss = 0
  do XYear = 1, WorkYearN
    WorkPredictor (XMem,XYear) = TimSeries (PredictorReg,XYear)
    
    if (WorkPredictor (XMem,XYear).EQ.MissVal) TotalMiss = TotalMiss + 1
  end do
  write (99,"(a23,i10)"), "predictand: total poss: ", WorkYearN 
  write (99,"(a23,i10)"), "predictor:  total miss: ", TotalMiss
  
  deallocate (TimRegNames,TimSeries)
  
  do
    print*, "  > PREDICTAND:"
    call LoadTim (WorkYearN,WorkADYear,TimRegN,TimRegNames,TimSeries)
    
    if (TimRegN.NE.WorkRegN) print*, "  > Unacceptable .tim file due to wrong no. of regions."
    if (TimRegN.EQ.WorkRegN) exit
  end do
  
  TotalMiss = 0
  do XReg = 1, WorkRegN
   do XYear = 1, WorkYearN
    WorkPredictand (XReg,XMem,XYear) = TimSeries (XReg,XYear)
    if (WorkPredictand (XReg,XMem,XYear).EQ.MissVal) TotalMiss = TotalMiss + 1
   end do
  end do
  write (99,"(a23,i10)"), "predictand: total poss: ", (WorkRegN*WorkYearN)  
  write (99,"(a23,i10)"), "predictand: total miss: ", TotalMiss  
  
  if (XMem.NE.WorkMemN) then
  	deallocate (TimRegNames,TimSeries)
  else
  	deallocate (TimSeries)			! we require TimRegNames for any saves to .tim
  end if
  
end do

end subroutine GetTimFiles

!*******************************************************************************
! calculate  'a' in y=ax

subroutine CalcAye

allocate (MemKey(WorkMemN),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure #####"
MemKey     = 0
SelectMemN = 0

if (WorkMemN.GT.1) then
  print*, "  > Enter the members to use in the calculation (99=exit): "
  do
	read (*,*,iostat=ReadStatus), SelectMem
	if (ReadStatus.LE.0.AND.SelectMem.GE.1.AND.SelectMem.LE.WorkMemN) &
		MemKey (SelectMem) = 1
	if (ReadStatus.LE.0.AND.SelectMem.EQ.99) exit
  end do
else
  MemKey (1) = 1
end if

do XMem = 1, WorkMemN
  if (MemKey(XMem).EQ.1) SelectMemN = SelectMemN + 1
end do

print*, "  > Calculate scaling eq using 30y slices ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), QEnsSlice
	if (ReadStatus.LE.0 .AND. QEnsSlice.GE.1.AND.QEnsSlice.LE.2) exit
end do
print*, "  > Calculating..."

if (AyeActive.EQ.1) then
  deallocate (WorkRegAye, stat=AllocStat)
  if (AllocStat.EQ.0) AyeActive = 0
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure #####"
end if

allocate (WorkActualExe  (SelectMemN,WorkYearN), &
	  WorkActualWye  (SelectMemN,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

allocate (WorkRegAye (WorkRegN), stat=AllocStat)
if (AllocStat.EQ.0) AyeActive = 1
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

XSelectMem = 0
do XMem = 1, WorkMemN
 if (MemKey(XMem).EQ.1) then
  XSelectMem = XSelectMem + 1
  do XYear = 1, WorkYearN
    WorkActualExe (XSelectMem,XYear) = WorkPredictor (XMem,XYear)
  end do
 end if
end do

if (QEnsSlice.EQ.2) then
  allocate (WorkOutExe(SelectMemN,WorkYearN), &
  	    WorkOutWye(SelectMemN,WorkYearN), &
  	    WorkInExe (WorkYearN),            &
  	    WorkInWye (WorkYearN),	      stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
  
  do XMem = 1, SelectMemN
   do XYear = 1, WorkYearN
    WorkOutExe (1,XYear) = WorkActualExe (XMem,XYear)
   end do
   WorkInExe = MissVal
   
   call Sliceify (1,WorkYearN,30,WorkMissAccept,WorkOutExe,WorkInExe)
   
   do XYear = 1, WorkYearN
    WorkActualExe (XMem,XYear) = WorkInExe (XYear)
   end do
  end do
end if
    
do XReg = 1, WorkRegN
  XSelectMem = 0
  do XMem = 1, WorkMemN
   if (MemKey(XMem).EQ.1) then
    XSelectMem = XSelectMem + 1
    do XYear = 1, WorkYearN
      WorkActualWye (XSelectMem,XYear) = WorkPredictand (XReg,XMem,XYear)
    end do
   end if
  end do
  
  if (QEnsSlice.EQ.2) then
   do XMem = 1, SelectMemN
    do XYear = 1, WorkYearN
     WorkOutWye (1,XYear) = WorkActualWye (XMem,XYear)
    end do
    WorkInWye = MissVal
   
    call Sliceify (1,WorkYearN,30,WorkMissAccept,WorkOutWye,WorkInWye)
   
    do XYear = 1, WorkYearN
     WorkActualWye (XMem,XYear) = WorkInWye (XYear)
    end do
   end do
  end if
  
  call LinearLSRZero (SelectMemN,WorkYearN,WorkActualExe,WorkActualWye,WorkAye)  
!  call LinearLSR (SelectMemN,WorkYearN,WorkActualExe,WorkActualWye,WorkAye,WorkBee,WorkPea)  

  WorkRegAye (XReg) = WorkAye
end do

deallocate (WorkActualExe,WorkActualWye,MemKey)
deallocate (WorkInExe,WorkInWye,WorkOutExe,WorkOutWye)

end subroutine CalcAye

!*******************************************************************************
! save region set of aye to .glo file

subroutine SaveAyeGlo

print*, "  > Enter the .glo file title (suggestions below): "
print*, trim(adjustl(WorkGridTitle)) // " : run(s) : period : " // trim(adjustl(WorkRegTitle)) &
  		// " : sc-eq 'a' : season "
do
	read (*,*,iostat=ReadStatus), WorkGloTitle
	if (ReadStatus.LE.0.AND.WorkGloTitle.NE."") exit
end do

call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGloTitle,WorkGridFilePath, &
	      WorkRegAye,WorkMapIDLReg)

end subroutine SaveAyeGlo

!*******************************************************************************
! load region set of aye from .glo file

subroutine LoadAyeGlo

if (AyeActive.EQ.1) then
  deallocate (WorkRegAye, stat=AllocStat)
  if (AllocStat.EQ.0) AyeActive = 0
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure #####"
end if

allocate (WorkRegAye (WorkRegN), stat=AllocStat)
if (AllocStat.EQ.0) AyeActive = 1
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

call LoadGlo (WorkLongN,WorkLatN,WorkRegN,WorkMapIDLReg,WorkRegAye)

end subroutine LoadAyeGlo

!*******************************************************************************
! make predictions using scaling eq on ensemble

subroutine ScalingPrediction

call PrepareToPredict
call SelectSaveRegions

XSelectMem = 0
do XMem = 1, WorkMemN		! load up predictor
 if (MemKey(XMem).EQ.1) then
  XSelectMem = XSelectMem + 1
  do XYear = 1, WorkYearN
    WorkActualExe (XSelectMem,XYear) = WorkPredictor (XMem,XYear)
  end do
 end if
end do

if (QEnsSlice.EQ.2) then
  allocate (WorkOutExe(SelectMemN,WorkYearN), &
  	    WorkOutWye(SelectMemN,WorkYearN), &
  	    WorkInExe (WorkYearN),            &
  	    WorkInWye (WorkYearN),	      stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
  
  do XMem = 1, SelectMemN
   write (99,"(a8,a20)"), "member: ", MemberName(XMem)

   do XYear = 1, WorkYearN
    WorkOutExe (1,XYear) = WorkActualExe (XMem,XYear)
   end do
   WorkInExe = MissVal
   
   call Sliceify (1,WorkYearN,30,WorkMissAccept,WorkOutExe,WorkInExe)
   
   TotalMiss = 0
   do XYear = 1, WorkYearN
    WorkActualExe (XMem,XYear) = WorkInExe (XYear)
    if (WorkActualExe(XMem,XYear).EQ.MissVal) TotalMiss = TotalMiss + 1
   end do
   write (99,"(a30,i10)"), "predictor: sliced: total poss: ", (WorkYearN*SelectMemN) 
   write (99,"(a30,i10)"), "predictor: sliced: total miss: ", TotalMiss
  end do
end if
    
TotalMiss = 0

do XReg = 1, WorkRegN		! loop by region
  XSelectMem = 0		! load up predictand 
  do XMem = 1, WorkMemN
   if (MemKey(XMem).EQ.1) then
    XSelectMem = XSelectMem + 1
    do XYear = 1, WorkYearN
      WorkActualWye (XSelectMem,XYear) = WorkPredictand (XReg,XMem,XYear)
    end do
   end if
  end do
  
  if (QEnsSlice.EQ.2) then
   do XMem = 1, SelectMemN
    do XYear = 1, WorkYearN
     WorkOutWye (1,XYear) = WorkActualWye (XMem,XYear)
    end do
    WorkInWye = MissVal
   
    call Sliceify (1,WorkYearN,30,WorkMissAccept,WorkOutWye,WorkInWye)
   
    do XYear = 1, WorkYearN
     WorkActualWye (XMem,XYear) = WorkInWye (XYear)
    end do
   end do
  end if
  
  WorkPredictWye  = MissVal	! initialise output arrays

  				! make predictions
  ParticAye = WorkRegAye(XReg)
  call PredictRun (SelectMemN,WorkYearN,WorkActualExe,WorkPredictWye,ParticAye)
  
  do XMem = 1, SelectMemN	! store pred and act
    do XYear = 1, WorkYearN
     if (WorkPredictWye(XMem,XYear).NE.MissVal.AND.WorkActualWye(XMem,XYear).NE.MissVal) then
      WorkAllPredict (XReg,XMem,XYear) = WorkPredictWye (XMem,XYear) 
      WorkAllActual  (XReg,XMem,XYear) = WorkActualWye  (XMem,XYear)
     else
      TotalMiss = TotalMiss + 1
     end if
    end do
  end do
  
  if (SaveErrStat.EQ.2) then	! if appropriate evaluate error of predictions
    				! the 3 at the front means calc Brier score
    WorkErrorMSE    = MissVal
    WorkInternalMSE = MissVal
    WorkFractionMSE = MissVal
    call PredictError (WorkMissAccept,3,SelectMemN,WorkYearN,WorkSliceLen,WorkActualWye,WorkPredictWye,&
    		       WorkErrorMSE,WorkInternalMSE,WorkFractionMSE)
    do XYear = 1, WorkYearN
      WorkAllFractionMSE (XReg,XYear) = WorkFractionMSE (XYear)
    end do
  end if
    				! if this region has been selected, then save a .lin
  if (SaveRegKey(XReg).EQ.1) call ActualSaveRegion    
end do
write (99,"(a34,i10)"), "predictand: predicted: total poss: ", (WorkRegN*WorkYearN*SelectMemN) 
write (99,"(a34,i10)"), "predictand: predicted: total miss: ", TotalMiss

if (SaveDiffStat.EQ.2) then
  print*, "  > OPPORTUNITIES to save diff stats to .tim and .glo and .lin"
  call DiffStatTimProcedure
  call DiffStatGloProcedure
  call DiffStatLinProcedure
  if (SelectMemN.GE.2) then
    print*, "  > 1. Calc intra-ensemble variability ? (1=no,2=yes)"
    do
	read (*,*,iostat=ReadStatus), QIntraVar
	if (ReadStatus.LE.0 .AND. QIntraVar.GE.1.AND.QIntraVar.LE.2) exit
    end do

    if (QIntraVar.EQ.2) then		
		call DiffStatNatTimProcedure
		call DiffStatNatGloProcedure
		call DiffStatNatLinProcedure
    end if
  end if
end if
if (SaveCorrStat.EQ.2) call CorrStatProcedure
if (SaveErrStat.EQ.2)  call ErrStatProcedure

deallocate (WorkActualExe,WorkActualWye,WorkPredictWye,SaveRegKey)
deallocate (WorkAllActual,WorkAllPredict)
deallocate (MemKey)
if (SaveErrStat.EQ.2)  deallocate (WorkErrorMSE,WorkInternalMSE,WorkFractionMSE, &
					WorkAllFractionMSE,WorkAllErrorStat)
if (SaveCorrStat.EQ.2) deallocate (WorkSpaceA,WorkSpaceB,WorkSpaceR,WorkTimeABR)

end subroutine ScalingPrediction

!*******************************************************************************
! prepare to predict using scaling equations

subroutine PrepareToPredict

allocate (MemKey(WorkMemN),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure #####"
MemKey     = 0
SelectMemN = 0

if (WorkMemN.GT.1) then
  print*, "  > Enter the members to estimate by scaling (99=exit): "
  do
	read (*,*,iostat=ReadStatus), SelectMem
	if (ReadStatus.LE.0.AND.SelectMem.GE.1.AND.SelectMem.LE.WorkMemN) &
		MemKey (SelectMem) = 1
	if (ReadStatus.LE.0.AND.SelectMem.EQ.99) exit
  end do
else
  MemKey (1) = 1
end if

do XMem = 1, WorkMemN
  if (MemKey(XMem).EQ.1) SelectMemN = SelectMemN + 1
end do

print*, "  > Predict using 30y slices ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), QEnsSlice
	if (ReadStatus.LE.0 .AND. QEnsSlice.GE.1.AND.QEnsSlice.LE.2) exit
end do

print*, "  > Select the stats to calc."

print*, "  > 1. Calc diffs between scaled and modelled? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), SaveDiffStat
	if (ReadStatus.LE.0.AND.SaveDiffStat.GE.1.AND.SaveDiffStat.LE.2) exit
end do

if (SaveDiffStat.EQ.2) then
  print*, "  > 1. Calc absolute diffs (=1) or RMSE (=2) ?"
  do
	read (*,*,iostat=ReadStatus), QGloRMSE
	if (ReadStatus.LE.0.AND.QGloRMSE.GE.1.AND.QGloRMSE.LE.2) exit
  end do
  
  if (SelectMemN.GE.2.AND.QGloRMSE.EQ.1) then
    print*, "  > 1. Calc using scatter of members (=1), or ens means (=2) ?"
    do
	read (*,*,iostat=ReadStatus), QDiffScat
	if (ReadStatus.LE.0.AND.QDiffScat.GE.1.AND.QDiffScat.LE.2) exit
    end do
  else
    QDiffScat = 1
  end if
end if

print*, "  > 2. Calc correlations between scaled and modelled? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), SaveCorrStat
	if (ReadStatus.LE.0.AND.SaveCorrStat.GE.1.AND.SaveCorrStat.LE.2) exit
end do

if (SelectMemN.GT.1) then
  print*, "  > 3. Calculate the mean Brier-based score? (1=no,2=yes)"
  do
	read (*,*,iostat=ReadStatus), SaveErrStat
	if (ReadStatus.LE.0.AND.SaveErrStat.GE.1.AND.SaveErrStat.LE.2) exit
  end do
else
  SaveErrStat = 1
end if

if (SaveErrStat.EQ.2) then
  print*, "  > 3. Select the slice length over which to calc means (for MSE): "
  do
	read (*,*,iostat=ReadStatus), WorkSliceLen
	if (ReadStatus.LE.0.AND.WorkSliceLen.GT.0) exit
  end do

  print*, "  > 3. Select period over which to average Brier-based score: "
  call SetUpErrorStat (WorkMissAccept,WorkADYear,WorkYearN,WorkSliceLen,&
  		       ErrStatMethod,ErrStatYear0,ErrStatYear1,&
  		       ErrStatStatThreshold,ErrStatMissThreshold)
end if

allocate (SaveRegKey (WorkRegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
SaveRegKey = 0

allocate (WorkAllPredict  (WorkRegN,SelectMemN,WorkYearN), &
	  WorkAllActual   (WorkRegN,SelectMemN,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
WorkAllActual  = MissVal
WorkAllPredict = MissVal

allocate (WorkActualExe  (SelectMemN,WorkYearN), &
	  WorkActualWye  (SelectMemN,WorkYearN), &
	  WorkPredictWye (SelectMemN,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
WorkActualExe  = MissVal
WorkActualWye  = MissVal
WorkPredictWye = MissVal

if (SaveErrStat.EQ.2) then
  allocate (WorkErrorMSE    (WorkYearN), &
	    WorkInternalMSE (WorkYearN), &
	    WorkFractionMSE (WorkYearN), &
	    WorkAllFractionMSE (WorkRegN,WorkYearN), &
	    WorkAllErrorStat   (WorkRegN),	     stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
  
  WorkErrorMSE       = MissVal
  WorkInternalMSE    = MissVal
  WorkFractionMSE    = MissVal
  WorkAllFractionMSE = MissVal
  WorkAllErrorStat   = MissVal
end if

if (SaveCorrStat.EQ.2) then
  allocate (WorkSpaceA (WorkRegN),    &
  	    WorkSpaceB (WorkRegN),    &
  	    WorkSpaceR (WorkRegN),    &
  	    WorkTimeABR (3,WorkYearN),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

  WorkSpaceA  = MissVal
  WorkSpaceB  = MissVal
  WorkSpaceR  = MissVal
  WorkTimeABR = MissVal
end if

end subroutine PrepareToPredict

!*******************************************************************************
! select regions to have a .lin file saved

subroutine SelectSaveRegions

print*, "  > 4. The results for any region may be saved to .lin file."
print*, "  > 4. Enter each region to be saved (0=list, -55=select all, -99=end):"
do
	read (*,*,iostat=ReadStatus), SaveRegSelect
	
	if (ReadStatus.LE.0) then
	  if (SaveRegSelect.EQ.0) then
	    do XReg = 1, WorkRegN
	      print "(I6,A40)", XReg, WorkRegNames(XReg)
	    end do
	  end if
	  
	  if (SaveRegSelect.GE.1.AND.SaveRegSelect.LE.WorkRegN) then
	    	SaveRegKey (SaveRegSelect) = 1
	  end if
	  
	  if (SaveRegSelect.EQ.-55) then
	    do XReg = 1, WorkRegN
	    	SaveRegKey (XReg) = 1
	    end do
	  end if
	end if
		
	if (ReadStatus.LE.0.AND.SaveRegSelect.LT.0) exit
end do

TotRegSelect = 0
do XReg = 1, WorkRegN
   if (SaveRegKey(XReg).EQ.1) TotRegSelect = TotRegSelect + 1
end do
print*, "  > Total regions to be individually saved: ", TotRegSelect

end subroutine SelectSaveRegions

!*******************************************************************************
! calc of diff stat

subroutine DiffStatTimProcedure

allocate (WorkOutTim (WorkRegN,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DiffStatTimProcedure: Allocation failure #####"
WorkOutTim = MissVal
    
print*, "  > 1. Save difference stats to .tim ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), QDiffTim
	if (ReadStatus.LE.0.AND.QDiffTim.GE.1.AND.QDiffTim.LE.2) exit
end do

if (QDiffTim.EQ.2) then
 MemThresh = SelectMemN * (100.0 - WorkMissAccept) / 100.0
 
 if      (QGloRMSE.EQ.1) then				! calc mean diff for each reg-year

  do XReg = 1, WorkRegN
    do XYear = 1, WorkYearN
      OpTotal = 0.0
      OpEn    = 0.0
      
      do XMem = 1, SelectMemN
        if (WorkAllPredict(XReg,XMem,XYear).NE.MissVal) then
          if (WorkAllActual(XReg,XMem,XYear).NE.MissVal) then
            OpTotal = OpTotal + WorkAllPredict(XReg,XMem,XYear) - WorkAllActual(XReg,XMem,XYear)
            OpEn    = OpEn    + 1.0
          end if
        end if
      end do
      
      if (OpEn.GT.MemThresh) then
      	WorkOutTim (XReg,XYear) = OpTotal / OpEn
      end if
    end do
  end do

 else if (QGloRMSE.EQ.2) then				! calc RMSE for each reg-year
 
  if (SelectMemN.GT.1) then
   print*, "  > Calculate for sample (=1) or population (=2) ?"
   do
	read (*,*,iostat=ReadStatus), QSampPop
	if (ReadStatus.LE.0.AND.QSampPop.GE.1.AND.QSampPop.LE.2) exit
   end do
  else
   QSampPop = 1
   print*, "  > We calc for sample, because n=1"
  end if
 
  do XReg = 1, WorkRegN
    do XYear = 1, WorkYearN
      OpTotal = 0.0
      OpEn    = 0.0
      
      do XMem = 1, SelectMemN
        if (WorkAllPredict(XReg,XMem,XYear).NE.MissVal) then
          if (WorkAllActual(XReg,XMem,XYear).NE.MissVal) then
            OpTotal = OpTotal + ((WorkAllPredict(XReg,XMem,XYear) - WorkAllActual(XReg,XMem,XYear)) ** 2)
            OpEn    = OpEn    + 1.0
          end if
        end if
      end do
      
      if (OpEn.GT.MemThresh) then
      	if (QSampPop.EQ.2) then
      	  if (OpEn.GT.1) then
            WorkOutTim (XReg,XYear) = OpTotal / (OpEn - 1)
            WorkOutTim (XReg,XYear) = sqrt ( WorkOutTim (XReg,XYear) )        	
      	  else
      	    WorkOutTim (XReg,XYear) = MissVal
      	  end if
      	else
            WorkOutTim (XReg,XYear) = OpTotal / OpEn
            WorkOutTim (XReg,XYear) = sqrt ( WorkOutTim (XReg,XYear) )        	
      	end if	
      end if
    end do
  end do

 end if
 
 print*, "  > Enter the .tim file title: "
 do
	read (*,*,iostat=ReadStatus), WorkGloTitle
	if (ReadStatus.LE.0.AND.WorkGloTitle.NE."") exit
 end do
    
 call SaveTim (WorkRegN, WorkYearN, WorkGloTitle, TimRegNames, WorkADYear, WorkOutTim)
 
end if

deallocate (WorkOutTim)

end subroutine DiffStatTimProcedure

!*******************************************************************************
! calc of diff stat

subroutine DiffStatGloProcedure

allocate (WorkOutReg (WorkRegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
    
do
  if (QGloRMSE.EQ.1) print*, "  > 1. Enter the year AD for which to calculate reg diffs (99=end): "
  if (QGloRMSE.EQ.2) print*, "  > 1. Enter the year AD for which to calculate reg RMSE (99=end): "
  do
	read (*,*,iostat=ReadStatus), SelectADYear
	if (ReadStatus.LE.0) exit
  end do
  
  SelectXYear = 0
  XYear       = 0
  do
    XYear = XYear + 1
    if (WorkADYear(XYear).EQ.SelectADYear) SelectXYear = XYear
    if (SelectXYear.GT.0.OR.XYear.EQ.WorkYearN) exit
  end do
  
  if (SelectXYear.GT.0) then
    WorkOutReg = MissVal
    
    do XReg = 1, WorkRegN
      RegTotal   = 0.0
      RegEn      = 0.0
      
      do XMem = 1, SelectMemN
        if (WorkAllPredict(XReg,XMem,SelectXYear).NE.MissVal &
        		.AND.WorkAllActual(XReg,XMem,SelectXYear).NE.MissVal) then
          if (QGloRMSE.EQ.1) then
            RegTotal = RegTotal + WorkAllPredict(XReg,XMem,SelectXYear) - WorkAllActual(XReg,XMem,SelectXYear)
            RegEn    = RegEn    + 1
          else if (QGloRMSE.EQ.2) then
            RegTotal =RegTotal+(WorkAllPredict(XReg,XMem,SelectXYear)-WorkAllActual(XReg,XMem,SelectXYear))**2
            RegEn    = RegEn    + 1
          end if
        end if
      end do
      
      if (RegEn.GE.1) then
      	if (QGloRMSE.EQ.1) then
      		WorkOutReg(XReg) = RegTotal / RegEn		! calcs abs diff
      	else if (QGloRMSE.EQ.2) then
      		WorkOutReg(XReg) = RegTotal / RegEn
      		WorkOutReg(XReg) = sqrt (WorkOutReg(XReg))	! calcs RMSE
      	end if
      end if
    end do
    
    print*, "  > Enter the .glo file title: "
    do
	read (*,*,iostat=ReadStatus), WorkGloTitle
	if (ReadStatus.LE.0.AND.WorkGloTitle.NE."") exit
    end do
    
    call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGloTitle,WorkGridFilePath, &
	          WorkOutReg,WorkMapIDLReg)
  
  else
    if (SelectADYear.NE.99) print*, "  > Unacceptable year AD."
  end if
  
  if (SelectADYear.EQ.99) exit
end do

deallocate (WorkOutReg)

end subroutine DiffStatGloProcedure

!*******************************************************************************
! calc of diff stat (global collection)

subroutine DiffStatLinProcedure

print*, "  > 1. Save glo RMSE to .lin ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), SelectDiffLin
	if (ReadStatus.LE.0.AND.SelectDiffLin.GE.1.AND.SelectDiffLin.LE.2) exit
end do

if (SelectDiffLin.EQ.2) then
  print*, "  > Calculating RMSE..."
  
  allocate (LinToBeSaved      (1,WorkYearN), &
    	    LinToBeSavedNames (1)          ,	stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"    
  LinToBeSaved          = MissVal
  LinToBeSavedNames (1) = "scaled v modelled : diff MSE across globe"
  
  if (QDiffScat.EQ.1) then
  	Threshold = real(WorkRegN) * real(SelectMemN) * (100.0 - WorkMissAccept) / 100.0
  else if (QDiffScat.EQ.2) then
  	Threshold = real(WorkRegN) * (100.0 - WorkMissAccept) / 100.0
  end if
  
  do XYear = 1, WorkYearN
    YearTotal = 0.0
    YearEn    = 0.0
    
    if (QDiffScat.EQ.1) then
      do XReg = 1, WorkRegN
       do XMem = 1, SelectMemN
        if (WorkAllPredict(XReg,XMem,XYear).NE.MissVal.AND.WorkAllActual(XReg,XMem,XYear).NE.MissVal) then
          YearTotal = YearTotal + (WorkAllPredict(XReg,XMem,XYear)-WorkAllActual(XReg,XMem,XYear)) ** 2
          YearEn    = YearEn    + 1.0
        end if
       end do
      end do
    else if (QDiffScat.EQ.2) then
      do XReg = 1, WorkRegN
       PredictTotal = 0.0
       ActualTotal  = 0.0
       MemEn        = 0.0
       
       do XMem = 1, SelectMemN
        if (WorkAllPredict(XReg,XMem,XYear).NE.MissVal.AND.WorkAllActual(XReg,XMem,XYear).NE.MissVal) then
         PredictTotal = PredictTotal + WorkAllPredict (XReg,XMem,XYear)
         ActualTotal  = ActualTotal  + WorkAllActual  (XReg,XMem,XYear)
         MemEn        = MemEn        + 1.0
        end if
       end do
       
       if (MemEn.GE.1) then
         YearTotal = YearTotal + ((PredictTotal/MemEn)-(ActualTotal/MemEn)) ** 2
         YearEn    = YearEn    + 1.0
       end if
      end do
    end if
    
    if (YearEn.GE.Threshold) then
    	LinToBeSaved(1,XYear) = YearTotal / YearEn
    	LinToBeSaved(1,XYear) = sqrt ( LinToBeSaved(1,XYear) )
    end if
  end do
  
  print*, "  > Enter the .lin file title: "
  do
	read (*,*,iostat=ReadStatus), SaveLinTitle
	if (ReadStatus.LE.0.AND.SaveLinTitle.NE."") exit
  end do

  call SaveLin (1,WorkYearN,SaveLinTitle,LinToBeSavedNames,WorkADYear,LinToBeSaved)
  
  deallocate (LinToBeSaved,LinToBeSavedNames)
end if

end subroutine DiffStatLinProcedure

!*******************************************************************************
! calc of .lin diff stat in ensemble

subroutine DiffStatNatLinProcedure

print*, "  > 1. Save RMSE of global nat var to .lin ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), SelectNatLin
	if (ReadStatus.LE.0.AND.SelectNatLin.GE.1.AND.SelectNatLin.LE.2) exit
end do

if (SelectNatLin.EQ.2) then
  print*, "  > Calculating RMSE..."
  
  allocate (LinToBeSaved      (1,WorkYearN), &
    	    LinToBeSavedNames (1)          ,	stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"    
  LinToBeSaved          = MissVal
  LinToBeSavedNames (1) = "modelled v modelled : nat MSE across globe"
  
  do XYear = 1, WorkYearN
    YearTotal = 0.0
    YearEn    = 0.0
    
    do XReg = 1, WorkRegN
      do XPredictor = 1, SelectMemN
        if (WorkAllActual(XReg,XPredictor,XYear).NE.MissVal) then
          do XPredictand = 1, SelectMemN
            if (XPredictor.NE.XPredictand.AND.WorkAllActual(XReg,XPredictand,XYear).NE.MissVal) then
              YearTotal = YearTotal + (WorkAllActual(XReg,XPredictor,XYear)- &
              				WorkAllActual(XReg,XPredictand,XYear)) ** 2
              YearEn    = YearEn    + 1.0
            end if
          end do
        end if
      end do
    end do

    if (YearEn.GE.2) then
    	LinToBeSaved(1,XYear) = YearTotal / YearEn
    	LinToBeSaved(1,XYear) = sqrt (LinToBeSaved(1,XYear))
    end if
  end do
  
  print*, "  > Enter the .lin file title: "
  do
	read (*,*,iostat=ReadStatus), SaveLinTitle
	if (ReadStatus.LE.0.AND.SaveLinTitle.NE."") exit
  end do

  call SaveLin (1,WorkYearN,SaveLinTitle,LinToBeSavedNames,WorkADYear,LinToBeSaved)
  
  deallocate (LinToBeSaved,LinToBeSavedNames)
end if

end subroutine DiffStatNatLinProcedure

!*******************************************************************************
! calc of .glo diff stat in ensemble

subroutine DiffStatNatGloProcedure

print*, "  > 1. Save RMSE of regional nat var to .glo ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), SelectNatLin
	if (ReadStatus.LE.0.AND.SelectNatLin.GE.1.AND.SelectNatLin.LE.2) exit
end do

if (SelectNatLin.EQ.2) then
  print*, "  > Calculating RMSE..."
  
  allocate (RegToBeSaved (WorkRegN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"    
  RegToBeSaved = MissVal
  
  do XReg = 1, WorkRegN
    RegTotal = 0.0
    RegEn    = 0.0
    
    do XYear = 1, WorkYearN
      do XPredictor = 1, SelectMemN
        if (WorkAllActual(XReg,XPredictor,XYear).NE.MissVal) then
          do XPredictand = 1, SelectMemN
            if (XPredictor.NE.XPredictand.AND.WorkAllActual(XReg,XPredictand,XYear).NE.MissVal) then
              RegTotal = RegTotal + (WorkAllActual(XReg,XPredictor,XYear)- &
              				WorkAllActual(XReg,XPredictand,XYear)) ** 2
              RegEn    = RegEn    + 1.0
            end if
          end do
        end if
      end do
    end do

    if (RegEn.GE.2) then
    	RegToBeSaved(XReg) = RegTotal / RegEn
    	RegToBeSaved(XReg) = sqrt (RegToBeSaved(XReg))
    end if      
  end do
  
  print*, "  > Enter the .glo file title: "
  do
	read (*,*,iostat=ReadStatus), WorkGloTitle
	if (ReadStatus.LE.0.AND.WorkGloTitle.NE."") exit
  end do

  call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGloTitle,WorkGridFilePath, &
	        RegToBeSaved,WorkMapIDLReg)
  
  deallocate (RegToBeSaved)
end if

end subroutine DiffStatNatGloProcedure

!*******************************************************************************
! calc of .tim diff stat in ensemble

subroutine DiffStatNatTimProcedure

print*, "  > 1. Save RMSE of nat var to .tim ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), SelectNatLin
	if (ReadStatus.LE.0.AND.SelectNatLin.GE.1.AND.SelectNatLin.LE.2) exit
end do

if (SelectNatLin.EQ.2) then
  print*, "  > Calculate for sample (=1) or population (=2) ?"
  do
	read (*,*,iostat=ReadStatus), QSampPop
	if (ReadStatus.LE.0.AND.QSampPop.GE.1.AND.QSampPop.LE.2) exit
  end do
 
  print*, "  > Calculating RMSE..."
  
  allocate (WorkOutTim (WorkRegN,WorkYearN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: DiffstatNatTimProcedure: Allocation failure #####"    
  WorkOutTim = MissVal
  
  do XReg = 1, WorkRegN
   do XYear = 1, WorkYearN
    OpTotal = 0.0
    OpEn    = 0.0
    
    do XPredictor = 1, SelectMemN
       if (WorkAllActual(XReg,XPredictor,XYear).NE.MissVal) then
          OpEn = OpEn + 1.0
          
          do XPredictand = 1, SelectMemN
            if (XPredictor.NE.XPredictand.AND.WorkAllActual(XReg,XPredictand,XYear).NE.MissVal) then
              OpTotal = OpTotal + (WorkAllActual(XReg,XPredictor,XYear)- &
              				WorkAllActual(XReg,XPredictand,XYear)) ** 2
            end if
          end do
       end if
    end do

    if (OpEn.GE.2) then
      	if (QSampPop.EQ.2) then
          WorkOutTim(XReg,XYear) = OpTotal / ((OpEn - 1) ** 2)
    	  WorkOutTim(XReg,XYear) = sqrt (WorkOutTim(XReg,XYear) )
      	else
          WorkOutTim(XReg,XYear) = OpTotal / (OpEn * (OpEn - 1))
    	  WorkOutTim(XReg,XYear) = sqrt (WorkOutTim(XReg,XYear) )
      	end if	
    end if
   end do
  end do
  
  print*, "  > Enter the .tim file title: "
  do
	read (*,*,iostat=ReadStatus), WorkGloTitle
	if (ReadStatus.LE.0.AND.WorkGloTitle.NE."") exit
  end do

  call SaveTim (WorkRegN, WorkYearN, WorkGloTitle, TimRegNames, WorkADYear, WorkOutTim)
   
  deallocate (WorkOutTim, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: DiffStatNatTimProcedure: Deallocation failure #####"    
end if

end subroutine DiffStatNatTimProcedure

!*******************************************************************************
! calc of correlation stats

subroutine CorrStatProcedure

allocate (WorkOutExe (SelectMemN,WorkYearN), &
	  WorkOutWye (SelectMemN,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

TotalMiss = 0
TotalMissAye = 0
TotalMissBee = 0
TotalMissPea = 0

do XReg = 1, WorkRegN		! calc spatially-varying corr
  do XMem = 1, SelectMemN
    do XYear = 1, WorkYearN
      WorkOutExe (XMem,XYear) = WorkAllActual  (XReg,XMem,XYear)
      WorkOutWye (XMem,XYear) = WorkAllPredict (XReg,XMem,XYear)
      
      if (WorkOutExe(XMem,XYear).EQ.MissVal.OR.WorkOutWye(XMem,XYear).EQ.MissVal) &
      		TotalMiss = TotalMiss + 1
    end do
  end do
  
  call LinearLSR (SelectMemN,WorkYearN,WorkOutExe,WorkOutWye,LSRAye,LSRBee,LSRPea)
  
  WorkSpaceA (XReg) = LSRAye
  WorkSpaceB (XReg) = LSRBee
  WorkSpaceR (XReg) = LSRPea
  
  if (LSRAye.EQ.MissVal) TotalMissAye = TotalMissAye + 1
  if (LSRBee.EQ.MissVal) TotalMissBee = TotalMissBee + 1
  if (LSRPea.EQ.MissVal) TotalMissPea = TotalMissPea + 1
end do

write (99,"(a33,i10)"), "corr stat data (glo): total poss: ", (WorkRegN*WorkYearN*SelectMemN) 
write (99,"(a33,i10)"), "corr stat data (glo): total miss: ", TotalMiss
write (99,"(a33,i10)"), "corr a,b,r     (glo): total poss: ", WorkRegN 
write (99,"(a33,3i10)"), "corr a,b,r     (glo): total miss: ", TotalMissAye, TotalMissBee, TotalMissPea 

deallocate (WorkOutExe,WorkOutWye)

allocate (WorkOutExe (SelectMemN,WorkRegN), &
	  WorkOutWye (SelectMemN,WorkRegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

do XYear = 1, WorkYearN		! calc temporally-varying corr
  do XMem = 1, SelectMemN
    do XReg = 1, WorkRegN
      WorkOutExe (XMem,XReg) = WorkAllActual  (XReg,XMem,XYear)
      WorkOutWye (XMem,XReg) = WorkAllPredict (XReg,XMem,XYear)
    end do
  end do
  
  call LinearLSR (SelectMemN,WorkRegN,WorkOutExe,WorkOutWye,LSRAye,LSRBee,LSRPea)
  
  WorkTimeABR (1,XYear) = LSRAye
  WorkTimeABR (2,XYear) = LSRBee
  WorkTimeABR (3,XYear) = LSRPea
end do

deallocate (WorkOutExe,WorkOutWye)

print*, "  > 2. Save glo stats to .lin ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), SaveTimeABR
	if (ReadStatus.LE.0.AND.SaveTimeABR.GE.1.AND.SaveTimeABR.LE.2) exit
end do

if (SaveTimeABR.EQ.2) then
  print*, "  > Enter the .lin file title: "
  do
	read (*,*,iostat=ReadStatus), SaveLinTitle
	if (ReadStatus.LE.0.AND.SaveLinTitle.NE."") exit
  end do
  
  allocate (LinToBeSavedNames (3), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
  LinToBeSavedNames(1) = "global LSR 'a' y=ax+b"
  LinToBeSavedNames(2) = "global LSR 'b' y=ax+b"
  LinToBeSavedNames(3) = "global Pearson coeff"
  
  call SaveLin (3,WorkYearN,SaveLinTitle,LinToBeSavedNames,WorkADYear,WorkTimeABR)
  
  deallocate (LinToBeSavedNames)
end if

print*, "  > 2. Auto-save reg LSR stats (y=ax+b) to 3 * .glo ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), SaveSpaceA
	if (ReadStatus.LE.0.AND.SaveSpaceA.GE.1.AND.SaveSpaceA.LE.2) exit
end do

if (SaveSpaceA.EQ.2) then
  print*, "  > Enter the generic .glo file title: "
  do
	read (*,*,iostat=ReadStatus), WorkGloTitle
	if (ReadStatus.LE.0.AND.WorkGloTitle.NE."") exit
  end do
  
  print*, "  > File to be saved includes: y=ax+b : a"
  call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGloTitle,WorkGridFilePath, &
	        WorkSpaceA,WorkMapIDLReg)
	        
  print*, "  > File to be saved includes: y=ax+b : b"
  call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGloTitle,WorkGridFilePath, &
	        WorkSpaceB,WorkMapIDLReg)
	        
  print*, "  > File to be saved includes: y=ax+b : r"
  call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGloTitle,WorkGridFilePath, &
	        WorkSpaceR,WorkMapIDLReg)	        
end if

end subroutine CorrStatProcedure

!*******************************************************************************
! calc of error stat

subroutine ErrStatProcedure

if (SaveErrStat.EQ.2) then
  call PredictErrorStat (WorkAllFractionMSE,WorkADYear,WorkRegN,WorkYearN,&
		         ErrStatMethod,ErrStatYear0,ErrStatYear1,&
		         ErrStatStatThreshold,ErrStatMissThreshold,WorkAllErrorStat)

  do XReg = 1, WorkRegN		! print out stat for chosen regions
    if (SaveRegKey(XReg).EQ.1) then    
      print*, "  > 3+4. Region no., mean Brier score: ", XReg, WorkAllErrorStat(XReg)
    end if
  end do

  print*, "  > 3. Enter the .glo file title: "
  do
	read (*,*,iostat=ReadStatus), WorkGloTitle
	if (ReadStatus.LE.0.AND.WorkGloTitle.NE."") exit
  end do
  
  TotalMiss = 0
  do XReg = 1, WorkRegN
      if (WorkAllErrorStat(XReg).EQ.MissVal) TotalMiss = TotalMiss + 1 
  end do
  write (9,"(a21,i10)"), "err stat: total miss: ", TotalMiss
  
  call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGloTitle,WorkGridFilePath, &
	        WorkAllErrorStat,WorkMapIDLReg)
end if

end subroutine ErrStatProcedure

!*******************************************************************************
! actual save region

subroutine ActualSaveRegion

print*, "  > 4. Save region: ", XReg, WorkRegNames (XReg)

if (SaveErrStat.EQ.2) then
  SaveLinN  = SelectMemN * 3 + 3
else
  SaveLinN  = SelectMemN * 3
end if

allocate (LinToBeSaved      (SaveLinN,WorkYearN), &
    	  LinToBeSavedNames (SaveLinN)          ,	stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
    
LinToBeSaved      = MissVal
LinToBeSavedNames = "blank"
    
if (SaveErrStat.EQ.2) then
      LinToBeSavedNames ((SelectMemN*3)+1) = "all runs: scaling MSE"
      LinToBeSavedNames ((SelectMemN*3)+2) = "all runs: control MSE"
      LinToBeSavedNames ((SelectMemN*3)+3) = "all runs: Brier score"

      do XYear = 1, WorkYearN
        LinToBeSaved (((SelectMemN*3)+1),XYear) = WorkErrorMSE    (XYear)
        LinToBeSaved (((SelectMemN*3)+2),XYear) = WorkInternalMSE (XYear)
        LinToBeSaved (((SelectMemN*3)+3),XYear) = WorkFractionMSE (XYear)
      end do
end if
    
XSelectMem = 0
do XMem = 1, WorkMemN
   if (MemKey(XMem).EQ.1) then
      XSelectMem = XSelectMem + 1
      
      LinToBeSavedNames ((SelectMemN*0)+XSelectMem) = trim(adjustl(MemberName(XMem))) // " " &
      			// "model predictor"
      LinToBeSavedNames ((SelectMemN*1)+XSelectMem) = trim(adjustl(MemberName(XMem))) // " " &
      			// "model predictand"
      LinToBeSavedNames ((SelectMemN*2)+XSelectMem) = trim(adjustl(MemberName(XMem))) // " " &
      			// "scaled predictand"

      do XYear = 1, WorkYearN
        LinToBeSaved (((SelectMemN*0)+XSelectMem),XYear) = WorkActualExe   (XSelectMem,XYear)
        LinToBeSaved (((SelectMemN*1)+XSelectMem),XYear) = WorkActualWye   (XSelectMem,XYear)
        LinToBeSaved (((SelectMemN*2)+XSelectMem),XYear) = WorkPredictWye  (XSelectMem,XYear)
      end do
   end if 
end do
    
GaussChoice = 0
print*, "  > Save as raw (=1) or after 30y Gaussian filtering (=2) ?"
do
	read (*,*,iostat=ReadStatus), GaussChoice
	if (ReadStatus.LE.0.AND.GaussChoice.GE.1.AND.GaussChoice.LE.2) exit
end do
    
if (GaussChoice.EQ.2) then
    allocate (TempUnSmoo  (WorkYearN), &
              TempLowSmoo (WorkYearN), &
      	      TempHighSmoo(WorkYearN), stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
      
    do XLine = 1, SaveLinN
        LinToBeSavedNames (XLine) = LinToBeSavedNames (XLine) // " 30GaussLow"
        
        TempUnSmoo   = MissVal
        TempLowSmoo  = MissVal
        TempHighSmoo = MissVal
        
        do XYear = 1, WorkYearN
          TempUnSmoo (XYear) = LinToBeSaved (XLine,XYear)
        end do
        
        call GaussSmooth (WorkYearN,30,TempUnSmoo,TempLowSmoo,TempHighSmoo)
        
        do XYear = 1, WorkYearN
          LinToBeSaved (XLine,XYear) = TempLowSmoo (XYear)
        end do        
    end do
      
    deallocate (TempUnSmoo,TempLowSmoo,TempHighSmoo)
end if
    
print*, "  > Enter the .lin file title (suggestions below): "
print*, " scale eq info : applic info "
print*, trim(adjustl(WorkGridTitle)) // " : run(s) : diagnostic : period : " &
    	    // " season : " // trim(adjustl(WorkRegTitle))    	    
do
	read (*,*,iostat=ReadStatus), SaveLinTitle
	if (ReadStatus.LE.0.AND.SaveLinTitle.NE."") exit
end do

call SaveLin (SaveLinN,WorkYearN,SaveLinTitle,LinToBeSavedNames,WorkADYear,LinToBeSaved)
    
deallocate (LinToBeSaved,LinToBeSavedNames,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure #####"
    
print*, "  > Saved."

end subroutine ActualSaveRegion

!*******************************************************************************
! calc nat var of current ensemble

subroutine PrepareCalcNatVar

allocate (MemKey(WorkMemN),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure #####"
MemKey     = 0
SelectMemN = 0

if (WorkMemN.GT.1) then
  print*, "  > Enter the members for which to calc nat var stats (99=exit): "
  do
	read (*,*,iostat=ReadStatus), SelectMem
	if (ReadStatus.LE.0.AND.SelectMem.GE.1.AND.SelectMem.LE.WorkMemN) &
		MemKey (SelectMem) = 1
	if (ReadStatus.LE.0.AND.SelectMem.EQ.99) exit
  end do
else
  MemKey (1) = 1
end if

do XMem = 1, WorkMemN
  if (MemKey(XMem).EQ.1) SelectMemN = SelectMemN + 1
end do

print*, "  > Predict using 30y slices ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), QEnsSlice
	if (ReadStatus.LE.0 .AND. QEnsSlice.GE.1.AND.QEnsSlice.LE.2) exit
end do

print*, "  > Preparing to calculate natural variability statistics..."

allocate (WorkAllPredict (WorkRegN,SelectMemN,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
WorkAllPredict = MissVal

XSelectMem = 0
do XMem = 1, WorkMemN		! load up predictor
 if (MemKey(XMem).EQ.1) then

  XSelectMem = XSelectMem + 1

  do XReg = 1, WorkRegN
   do XYear = 1, WorkYearN
    WorkAllPredict (XReg,XSelectMem,XYear) = WorkPredictand (XReg,XMem,XYear)
   end do
  end do
  
 end if
end do

if (QEnsSlice.EQ.2) then
 allocate (WorkOutExe(SelectMemN,WorkYearN), &
  	   WorkInExe (WorkYearN),            stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
 
 do XReg = 1, WorkRegN 
  do XMem = 1, SelectMemN

   do XYear = 1, WorkYearN
    WorkOutExe (1,XYear) = WorkAllPredict (XReg,XMem,XYear)
   end do
   WorkInExe = MissVal
   
   call Sliceify (1,WorkYearN,30,WorkMissAccept,WorkOutExe,WorkInExe)
   
   do XYear = 1, WorkYearN
    WorkAllPredict (XReg,XMem,XYear) = WorkInExe (XYear)
   end do
  
  end do
 end do

 deallocate (WorkOutExe,WorkInExe)
end if
    
allocate (WorkSpaceABR (6,WorkRegN),    &
          WorkSpaceA   (WorkRegN), 	&
  	  WorkTimeABR  (6,WorkYearN),   stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

WorkSpaceABR = MissVal
WorkTimeABR  = MissVal

end subroutine PrepareCalcNatVar

!*******************************************************************************
! saves nat var stats

subroutine SaveNatVar

print*, "  > Save globally averaged stats to .lin ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), SaveTimeABR
	if (ReadStatus.LE.0.AND.SaveTimeABR.GE.1.AND.SaveTimeABR.LE.2) exit
end do

if (SaveTimeABR.EQ.2) then
  print*, "  > Enter the .lin file title: "
  do
	read (*,*,iostat=ReadStatus), SaveLinTitle
	if (ReadStatus.LE.0.AND.SaveLinTitle.NE."") exit
  end do
  
  allocate (LinToBeSavedNames (6), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
  LinToBeSavedNames(1) = "global LSR 'a' y=ax+b : min"
  LinToBeSavedNames(2) = "global LSR 'a' y=ax+b : max"
  LinToBeSavedNames(3) = "global LSR 'b' y=ax+b : min"
  LinToBeSavedNames(4) = "global LSR 'b' y=ax+b : min"
  LinToBeSavedNames(5) = "global Pearson coeff  : min"
  LinToBeSavedNames(6) = "global Pearson coeff  : max"
  
  call SaveLin (6,WorkYearN,SaveLinTitle,LinToBeSavedNames,WorkADYear,WorkTimeABR)
  
  deallocate (LinToBeSavedNames)
end if


print*, "  > Auto-save regionally av LSR stats (y=ax+b) to 6 * .glo ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), SaveSpaceA
	if (ReadStatus.LE.0.AND.SaveSpaceA.GE.1.AND.SaveSpaceA.LE.2) exit
end do

if (SaveSpaceA.EQ.2) then
  print*, "  > Enter the generic .glo file title: "
  do
	read (*,*,iostat=ReadStatus), WorkGloTitle
	if (ReadStatus.LE.0.AND.WorkGloTitle.NE."") exit
  end do
  
  do XGloFile = 1, 6
    print*, "  > File to be saved includes: "
    if (XGloFile.EQ.1.OR.XGloFile.EQ.2) print*, "  > y=ax+b : a"
    if (XGloFile.EQ.3.OR.XGloFile.EQ.4) print*, "  > y=ax+b : b"
    if (XGloFile.EQ.5.OR.XGloFile.EQ.6) print*, "  > y=ax+b : r"
    if (XGloFile.EQ.1.OR.XGloFile.EQ.3.OR.XGloFile.EQ.5) print*, "  > ensemble min values"
    if (XGloFile.EQ.2.OR.XGloFile.EQ.4.OR.XGloFile.EQ.6) print*, "  > ensemble max values"
      
    do XReg = 1, WorkRegN
      WorkSpaceA (XReg) = WorkSpaceABR (XGloFile,XReg)
    end do
  
    call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGloTitle,WorkGridFilePath, &
	          WorkSpaceA,WorkMapIDLReg)
  end do
end if

end subroutine SaveNatVar

!*******************************************************************************
! calcs nat var stats

subroutine CalcNatVar

call PrepareCalcNatVar

print*, "  > Calculating spatially-varying correlations..."

allocate (WorkOutExe (1,WorkYearN), &
	  WorkOutWye (1,WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

do XReg = 1, WorkRegN		! calc spatially-varying corr
  do XPredictor = 1, SelectMemN
   do XPredictand = 1, SelectMemN
    if (XPredictor.NE.XPredictand) then
     
     do XYear = 1, WorkYearN
      WorkOutExe (1,XYear) = WorkAllPredict (XReg,XPredictor, XYear)
      WorkOutWye (1,XYear) = WorkAllPredict (XReg,XPredictand,XYear)
     end do

     call LinearLSR (1,WorkYearN,WorkOutExe,WorkOutWye,LSRAye,LSRBee,LSRPea)
     
     if (LSRAye.NE.MissVal) then
       if (WorkSpaceABR(1,XReg).EQ.MissVal) WorkSpaceABR(1,XReg) = LSRAye
       if (WorkSpaceABR(2,XReg).EQ.MissVal) WorkSpaceABR(2,XReg) = LSRAye
       if (WorkSpaceABR(1,XReg).GT.LSRAye)  WorkSpaceABR(1,XReg) = LSRAye
       if (WorkSpaceABR(2,XReg).LT.LSRAye)  WorkSpaceABR(2,XReg) = LSRAye
     end if
     
     if (LSRBee.NE.MissVal) then
       if (WorkSpaceABR(3,XReg).EQ.MissVal) WorkSpaceABR(3,XReg) = LSRBee
       if (WorkSpaceABR(4,XReg).EQ.MissVal) WorkSpaceABR(4,XReg) = LSRBee
       if (WorkSpaceABR(3,XReg).GT.LSRBee)  WorkSpaceABR(3,XReg) = LSRBee
       if (WorkSpaceABR(4,XReg).LT.LSRBee)  WorkSpaceABR(4,XReg) = LSRBee
     end if
     
     if (LSRPea.NE.MissVal) then
       if (WorkSpaceABR(5,XReg).EQ.MissVal) WorkSpaceABR(5,XReg) = LSRPea
       if (WorkSpaceABR(6,XReg).EQ.MissVal) WorkSpaceABR(6,XReg) = LSRPea
       if (WorkSpaceABR(5,XReg).GT.LSRPea)  WorkSpaceABR(5,XReg) = LSRPea
       if (WorkSpaceABR(6,XReg).LT.LSRPea)  WorkSpaceABR(6,XReg) = LSRPea
     end if
     
    end if
   end do
  end do
end do

deallocate (WorkOutExe,WorkOutWye)

print*, "  > Calculating temporally-varying correlations..."

allocate (WorkOutExe (1,WorkRegN), &
	  WorkOutWye (1,WorkRegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

do XYear = 1, WorkYearN		! calc temporally-varying corr
  do XPredictor = 1, SelectMemN
   do XPredictand = 1, SelectMemN
    if (XPredictor.NE.XPredictand) then

     do XReg = 1, WorkRegN
      WorkOutExe (1,XReg) = WorkAllPredict (XReg,XPredictor, XYear)
      WorkOutWye (1,XReg) = WorkAllPredict (XReg,XPredictand,XYear)
     end do

     call LinearLSR (1,WorkYearN,WorkOutExe,WorkOutWye,LSRAye,LSRBee,LSRPea)
     
     if (LSRAye.NE.MissVal) then
       if (WorkTimeABR(1,XYear).EQ.MissVal) WorkTimeABR(1,XYear) = LSRAye
       if (WorkTimeABR(2,XYear).EQ.MissVal) WorkTimeABR(2,XYear) = LSRAye
       if (WorkTimeABR(1,XYear).GT.LSRAye)  WorkTimeABR(1,XYear) = LSRAye
       if (WorkTimeABR(2,XYear).LT.LSRAye)  WorkTimeABR(2,XYear) = LSRAye
     end if
     
     if (LSRBee.NE.MissVal) then
       if (WorkTimeABR(3,XYear).EQ.MissVal) WorkTimeABR(3,XYear) = LSRBee
       if (WorkTimeABR(4,XYear).EQ.MissVal) WorkTimeABR(4,XYear) = LSRBee
       if (WorkTimeABR(3,XYear).GT.LSRBee)  WorkTimeABR(3,XYear) = LSRBee
       if (WorkTimeABR(4,XYear).LT.LSRBee)  WorkTimeABR(4,XYear) = LSRBee
     end if
     
     if (LSRPea.NE.MissVal) then
       if (WorkTimeABR(5,XYear).EQ.MissVal) WorkTimeABR(5,XYear) = LSRPea
       if (WorkTimeABR(6,XYear).EQ.MissVal) WorkTimeABR(6,XYear) = LSRPea
       if (WorkTimeABR(5,XYear).GT.LSRPea)  WorkTimeABR(5,XYear) = LSRPea
       if (WorkTimeABR(6,XYear).LT.LSRPea)  WorkTimeABR(6,XYear) = LSRPea
     end if
     
    end if
   end do
  end do
end do

deallocate (WorkOutExe,WorkOutWye,WorkAllPredict)

call SaveNatVar

deallocate (MemKey,WorkSpaceABR,WorkTimeABR,WorkSpaceA)

end subroutine CalcNatVar

!*******************************************************************************
! close operations down

subroutine WindDown

if (AyeActive.EQ.1) then
  deallocate (WorkRegAye, stat=AllocStat)
  if (AllocStat.EQ.0) AyeActive = 0
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure #####"
end if

if (PredictActive.EQ.1) then
  deallocate (WorkPredictor,WorkPredictand,stat=AllocStat)
  if (AllocStat.EQ.0) PredictActive = 0
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure #####"
end if

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

print*

end subroutine WindDown

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

end program ScaleOps