! samplingerror2.f90
! f90 main program written on 16.02.99 by Tim Mitchell
! last modification on 17.07.00
! f90 -o samplingerror2 initialmod.f90 runselectmod.f90 extractmod.f90 scalemod.f90 savemod.f90 loadmod.f90
!        transformmod.f90 sortmod.f90 basicfun.f90 samplingerror2.f90

program SamplingError2

use InitialMod
use RunSelectMod
use ExtractMod
use ScaleMod
use SaveMod
use LoadMod
use TransformMod
use SortMod
use BasicFun

implicit none

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

real, pointer, dimension (:)	:: WorkGloAnaSlice, AllRMSE, PerGlo, ValidGlo
real, pointer, dimension (:)	:: AllAye, AllBee, OpOut, LowerBound, UpperBound, Weights
real, pointer, dimension (:,:)	:: WorkGotDecade, WorkGotFull, OpProc
real, pointer, dimension (:,:)	:: MyYear, MyAnom

integer, pointer, dimension (:) 	:: WorkADYear, WorkMapRawReg, WorkRegSizes
integer, pointer, dimension (:) 	:: WorkDecYearN, WorkDecGetYear0, WorkDecGetYear1
integer, pointer, dimension (:)		:: WorkDecStyleThis, WorkDecStyleNext
integer, pointer, dimension (:,:) 	:: WorkMapIDLRaw, WorkMapIDLReg
integer, pointer, dimension (:,:) 	:: WorkDecBlockKey

character (len=20), pointer, dimension (:) 	:: WorkRegNames
character (len=80), pointer, dimension (:) 	:: WorkDecPathThis, WorkDecPathNext

real, allocatable, dimension (:,:)	:: WorkBasePerMeans, WorkAnomPerMeans

integer, allocatable, dimension (:)	:: BaseSelect, PertSelect, UsePer

integer, dimension (1)	:: Seed

real, parameter :: MissVal = -999.0

real :: MissAccept, TotalValid, PerCentValid
real :: PerTot, PerEn, PerSqTot, RegTot, RegEn, RegSqTot
real :: PertThreshold, BaseThreshold, RMSEThreshold, RegThresh, PerThresh, YearThresh, RegPerThresh
real :: MemThresh, TrialThresh
real :: EnRMSE, TotRMSE
real :: MyAye, MyBee, MyPea
real :: ValidSqSum, RegSignal
real :: NoiseJFB, NoiseJFBPop
real :: RandomReal, CalcReal, SigLevel

integer :: WorkMonth0,WorkMonth1,WorkMonthN,WorkYearN,WorkDecN
integer :: WorkGrid,WorkLongN,WorkLatN,WorkDataN,WorkRegN
integer :: AllocStat, ReadStatus
integer :: XDec, XReg, XYear, XMonth, XPer, XLong, XTrial, XMem
integer :: Year0, Year1, WorkPerN
integer :: QAnother, QOrigBase, QAnomType, QGetTrend, QSaveLSR, QSampPop, QGetWeights
integer :: PerLen, GapLen, PerYear0, PerYear1, OrigLen, OrigBaseLen, PerYear, ThisPer
integer :: ValidReg, ValidYear, ValidCheck
integer :: ThreshExceeded
integer :: PerValid, RegValid, RegPerValid
integer :: WorkTrialN, RandomMem, CalcInt, TailN, RegMiss, TotMiss
integer :: EnsPerN, EnsPerYearN, EnsMemN

character (len=10) :: WorkGridTitle
character (len=40) :: WorkRegTitle
character (len=80) :: WorkGridFilePath, WorkDecTitle, SaveTitle

!*******************************************************************************
! main loop

MissAccept = 10.0

call Preliminaries
call GrabRun

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

do
  print*, "  >  0. Change % acceptable missing (currently: ", MissAccept
  print*, "  >  1. Detrend the raw"
  print*, "  >  2. Calc sd of anomaly means"
  print*, "  >  3. Calc pattern noise following JFB"
  print*, "  >  4. Calc sd of sd (interdec?)"
  print*, "  >  5. Calc interann var for ens of 30y slices: save mean+sd of set"
  print*, "  > 10. Check validity of file data and output into .glo"
  print*, "  > 99. Exit"
  
  do
	read (*,*,iostat=ReadStatus), QAnother
	if (ReadStatus.LE.0 .AND. QAnother.GE.0 .AND. QAnother.LE.99) exit
  end do
  
  if      (QAnother.EQ.0) then
  	call SetMissAccept (MissAccept)
  else if (QAnother.EQ.1) then
  	call DeTrendWork
  else if (QAnother.EQ.2) then
  	call ControlOne
  else if (QAnother.EQ.3) then
  	call ControlOne
  else if (QAnother.EQ.4) then
  	call ControlOne
  else if (QAnother.EQ.5) then
  	call InterAnnVar
  else if (QAnother.EQ.10) then
  	call CheckValid
  end if
    
  if (QAnother.EQ.99) exit
end do

print*
deallocate (WorkGotFull)
deallocate (WorkMapIDLRaw,WorkMapRawReg,WorkMapIDLReg,WorkRegSizes,WorkRegNames)
deallocate (PertSelect,BaseSelect)
close (99)

contains

!*******************************************************************************
! check validity of basic data array

subroutine CheckValid

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

YearThresh = WorkYearN * (100.0 - MissAccept) / 100.0
ValidReg   = 0
ValidCheck = 0

do XReg = 1, WorkRegN			! check contents of main data array
  ValidYear  = 0
  ValidSqSum = 0.0

  do XYear = 1, WorkYearN
    if (WorkGotFull(XReg,XYear).NE.MissVal) then
    	ValidYear  = ValidYear  + 1
    	ValidSqSum = ValidSqSum + WorkGotFull(XReg,XYear) ** 2
    end if
  end do
  
  ValidGlo (XReg) = ValidYear
  
  if (ValidYear .GE.YearThresh) then
  	ValidReg = ValidReg + 1
  	if (ValidSqSum.EQ.0.0) then
  		ValidCheck = ValidCheck + 1
  		write (99,"(a23,i6)"), "##### all zeros in reg: ", XReg
  	end if
  end if

  write (99,"(2i6,e12.5)"), XReg, ValidYear, ValidSqSum
end do

if (ValidCheck.GE.1) then
  print*, "  > ##### Suspect data: all zeros in regions, total: ", ValidCheck
  print*, "  > ##### See logsamerr.txt for details #####"
end if
print*, "  > We have valid regions totalling: ", ValidReg

call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,ValidGlo,WorkMapIDLReg)

end subroutine CheckValid

!*******************************************************************************
! controls calculation of options 2 and 3 and 4

subroutine ControlOne

print*, "  > Calc for absolute (=1) or percentage (=2) anomalies ?"
do
	read (*,*,iostat=ReadStatus), QAnomType
	if (ReadStatus.LE.0 .AND. QAnomType.GE.1 .AND. QAnomType.LE.2) exit
end do

if (WorkYearN.GT.240) then
     print*, "  > Is the original base the 240y control (=1) or a perturbed period (=2) ?"
     do
	read (*,*,iostat=ReadStatus), QOrigBase
	if (ReadStatus.LE.0 .AND. QOrigBase.GE.1 .AND. QOrigBase.LE.2) exit
     end do
else
     print*, "  > We assume that the original base is a perturbed period."
     QOrigBase = 2
end if
    
if (QOrigBase.EQ.1) then
      call Specification1
      call IdentifyPer1
      call CalcBaseMeans1
else if (QOrigBase.EQ.2) then
      call Specification2
      call IdentifyPer2
      call CalcBaseMeans2
end if

call CalcAnomMeans        

if      (QAnother.EQ.2) then
  call CalcStDev
else if (QAnother.EQ.3) then
  call GetWeights
  call CalcPerGlo
  call CalcJFBGlo
else if (QAnother.EQ.4) then
  call CalcStDevStDev
end if

deallocate (WorkAnomPerMeans)

end subroutine ControlOne

!*******************************************************************************
! loading prec preliminaries

subroutine Preliminaries

print*
print*, "  > ***** SamplingError *****"
print*, "  > Calculates sampling error from control data."
print*

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

allocate ( WorkGotDecade   (WorkRegN, 10), &
	   WorkGotFull     (WorkRegN, WorkYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"

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

print*, "  > Select the control run."

call RunSelect    (WorkGrid,WorkMonth0,WorkMonth1,WorkYearN,WorkDecN,WorkADYear,&
		   WorkDecTitle,WorkDecStyleThis,WorkDecStyleNext,&
		   WorkDecPathThis,WorkDecPathNext,WorkDecYearN,WorkDecGetYear0,WorkDecGetYear1)

end subroutine Preliminaries

!*******************************************************************************
! extract data into WorkGotFull

subroutine GrabRun

WorkGotFull = 0.0
print*, "  > Operating on decade starting in: "

do XDec = 1, WorkDecN
  Year0 = (XDec-1)*10 + 1
  Year1 = Year0 + 9
  
  WorkGotDecade = 0.0
  print*, WorkADYear(Year0)
  
  call BlockKey    (WorkDecYearN(XDec),WorkMonthN,WorkDecGetYear0(XDec),WorkDecGetYear1(XDec),&
  		    WorkMonth0,WorkMonth1,WorkDecPathThis(XDec),WorkDecPathNext(XDec),&
  		    WorkDecStyleThis(XDec),WorkDecStyleNext(XDec),WorkDecBlockKey)
  
  call ExtractFile (WorkLongN,WorkLatN,WorkDataN,WorkDecPathThis(XDec),WorkDecPathNext(XDec),&
  		    WorkDecStyleThis(XDec),WorkDecStyleNext(XDec),WorkRegN,WorkMonthN,&
  		    WorkDecYearN(XDec),WorkDecGetYear0(XDec),&
  		    WorkMapRawReg,WorkRegSizes,WorkDecBlockKey,WorkGotDecade)
  
  do XReg = 1, WorkRegN
    WorkGotFull (XReg,Year0:Year1) = WorkGotDecade (XReg,1:10)
  end do
end do

deallocate (WorkGotDecade,WorkDecYearN,WorkDecGetYear0,WorkDecGetYear1)
deallocate (WorkDecPathThis,WorkDecPathNext,WorkDecStyleThis,WorkDecStyleNext,WorkDecBlockKey)

if (WorkLongN.EQ.96.AND.WorkLatN.EQ.73.AND.WorkRegN.EQ.7008) then
  print*, "  > Converting regperbox into rpb..."
  
  do XLong = 1, WorkLongN
    XReg = WorkMapIDLReg (XLong,1)
    WorkGotFull (XReg,1:WorkYearN) = MissVal
    XReg = WorkMapIDLReg (XLong,WorkLatN)
    WorkGotFull (XReg,1:WorkYearN) = MissVal
  end do
end if

end subroutine GrabRun

!*******************************************************************************
! detrend raw data

subroutine DeTrendWork

call DeTrendTim (WorkRegN,WorkYearN,MissAccept,1,WorkADYear,WorkGotFull,AllAye,AllBee)

print*, "  > Save the values of 'a' and 'b' to .glo ? (1=no,2=yes)"
do
  read (*,*,iostat=ReadStatus), QSaveLSR
  if (ReadStatus.LE.0 .AND. QSaveLSR.GE.1 .AND. QSaveLSR.LE.2) exit
end do

if (QSaveLSR.EQ.2) then
  do XReg = 1, WorkRegN
    if (AllAye(XReg).NE.MissVal) AllAye(XReg) = AllAye(XReg) * 100.0
  end do
  
  print*, "  > Save 'a' as 100*trend-per-yearAD: "
  call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,AllAye,WorkMapIDLReg)

  print*, "  > Save 'b' (intercept): "
  call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,AllBee,WorkMapIDLReg)
end if

deallocate (AllAye,AllBee)

end subroutine DeTrendWork

!*******************************************************************************
! specify criteria for calc of sampling error

subroutine Specification1

print*, "  > Enter the period length for which to calc sampling error: "
do
	read (*,*,iostat=ReadStatus), PerLen
	if (ReadStatus.LE.0 .AND. PerLen.GE.1) exit
end do

print*, "  > Enter the gap length to leave between periods: "
do
	read (*,*,iostat=ReadStatus), GapLen
	if (ReadStatus.LE.0 .AND. GapLen.GE.0) exit
end do  

write (99,*), "specified criteria 1"

end subroutine Specification1

!*******************************************************************************
! specify criteria for calc of sampling error

subroutine Specification2

print*, "  > Enter the length used (base+anom) of the original: "
do
	read (*,*,iostat=ReadStatus), OrigLen
	if (ReadStatus.LE.0 .AND. OrigLen.GE.1) exit
end do

print*, "  > Enter the length of the original base: "
do
	read (*,*,iostat=ReadStatus), OrigBaseLen
	if (ReadStatus.LE.0 .AND. OrigBaseLen.GE.1 .AND. OrigBaseLen.LT.OrigLen) exit
end do  
  
print*, "  > Enter the first, last index yrs of the original to calc: (min,max poss: "
print "(2i6)", (OrigBaseLen+1), OrigLen
do
	read (*,*,iostat=ReadStatus), PerYear0, PerYear1
	if (ReadStatus.LE.0 .AND. PerYear1.GE.PerYear0 .AND. &
		PerYear0.GE.(OrigBaseLen+1) .AND. PerYear1.LE.OrigLen) exit
end do    

write (99,*), "specified criteria 2"

end subroutine Specification2

!*******************************************************************************
! identify periods

subroutine IdentifyPer1

BaseSelect = 0
PertSelect = 0

do XYear = 1, 240
    BaseSelect(XYear) = 1
end do

WorkPerN = floor ( (real(WorkYearN)-240.0) / (real(PerLen)+real(GapLen)) )
print*, "  > Maximum periods available: ", WorkPerN

do XPer = 1, WorkPerN
    Year0 = 241 + (XPer-1)*(PerLen+GapLen)
    do XYear = Year0, (Year0+PerLen-1)
      PertSelect (XYear) = XPer
    end do
end do

write (99,*), "identified periods 1"

end subroutine IdentifyPer1

!*******************************************************************************
! identify periods

subroutine IdentifyPer2

BaseSelect = 0
PertSelect = 0

WorkPerN = floor ( real(WorkYearN) / real(OrigLen) )
print*, "  > Maximum periods available: ", WorkPerN

do XPer = 1, WorkPerN
    Year0 = 1 + (XPer-1)*OrigLen
    do XYear = Year0, (Year0+OrigBaseLen-1)
      BaseSelect(XYear) = XPer
    end do
    
    Year0 = Year0 + PerYear0 - 1
    Year1 = Year0 + PerYear1 - PerYear0
    do XYear = Year0, Year1
      PertSelect(XYear) = XPer
    end do
end do

write (99,*), "identified periods 2"

end subroutine IdentifyPer2

!*******************************************************************************
! identify periods

subroutine IdentifyPer3

PertSelect = 0

WorkPerN = floor ( real(WorkYearN) / (real(PerLen)+real(GapLen)) )
print*, "  > Maximum periods available: ", WorkPerN

do XPer = 1, WorkPerN
    Year0 = (XPer-1)*(PerLen+GapLen)
    do XYear = Year0, (Year0+PerLen-1)
      PertSelect (XYear) = XPer
    end do
end do

end subroutine IdentifyPer3

!*******************************************************************************
! calc base period means

subroutine CalcBaseMeans1

allocate ( WorkBasePerMeans (WorkRegN,WorkPerN), &
 	   WorkAnomPerMeans (WorkRegN,WorkPerN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure #####"
WorkBasePerMeans = MissVal
WorkAnomPerMeans = MissVal

TotalValid = 0.0    
BaseThreshold = 240.0 *(100.0-MissAccept)/100.0
PertThreshold = PerLen*(100.0-MissAccept)/100.0

write (99,"(A16,F10.4)"), "base threshold: ", BaseThreshold
write (99,"(A16,F10.4)"), "pert threshold: ", PertThreshold
    
do XReg = 1, WorkRegN
    PerTot = 0.0
    PerEn  = 0.0
    
    do XYear = 1, WorkYearN
      if (BaseSelect(XYear).EQ.1) then
        if (WorkGotFull(XReg,XYear).NE.MissVal) then
          PerTot = PerTot + WorkGotFull(XReg,XYear)
          PerEn  = PerEn  + 1.0
        end if
      end if
    end do
    
!    write (99,*), "Region, Total-in-base, Years-in-base"
!    write (99,"(i6,2f10.4)"), XReg, PerTot, PerEn
    
    if (PerEn.GE.BaseThreshold.AND.PerTot.NE.0) then
      do XPer = 1, WorkPerN
        WorkBasePerMeans(XReg,XPer) = PerTot / PerEn
        TotalValid = TotalValid + 1.0
      end do
    end if
end do

PerCentValid = 100.0 * TotalValid / real(WorkRegN*WorkPerN)
print*, "  > Percentage period-region bases valid: ", PerCentValid
    
end subroutine CalcBaseMeans1

!*******************************************************************************
! calc base period means

subroutine CalcBaseMeans2

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

TotalValid = 0.0    
WorkBasePerMeans = MissVal
WorkAnomPerMeans = MissVal

BaseThreshold = OrigBaseLen*(100.0-MissAccept)/100.0
PertThreshold = (PerYear1-PerYear0+1.0)*(100.0-MissAccept)/100.0
    
write (99,"(A16,F10.4)"), "base threshold: ", BaseThreshold
write (99,"(A16,F10.4)"), "pert threshold: ", PertThreshold
    
do XReg = 1, WorkRegN
   do XPer = 1, WorkPerN
    PerTot = 0.0
    PerEn  = 0.0
    
    do XYear = 1, WorkYearN
      if (BaseSelect(XYear).EQ.XPer) then
        if (WorkGotFull(XReg,XYear).NE.MissVal) then
          PerTot = PerTot + WorkGotFull(XReg,XYear)
          PerEn  = PerEn  + 1.0
        end if
      end if
    end do
    
!    write (99,"(2i6,2f10.4)"), XReg, XPer, PerTot, PerEn
    
    if (PerEn.GE.BaseThreshold.AND.PerTot.NE.0) then
      WorkBasePerMeans(XReg,XPer) = PerTot / PerEn
      TotalValid = TotalValid + 1.0
    end if
   end do 
end do
      
PerCentValid = 100.0 * TotalValid / real(WorkRegN*WorkPerN)
print*, "  > Percentage period-region bases valid: ", PerCentValid
    
end subroutine CalcBaseMeans2

!*******************************************************************************
! calc anom period means

subroutine CalcAnomMeans

TotalValid = 0.0

!write (99,*), "in groups of three lines to one region and period"

do XReg = 1, WorkRegN
 do XPer = 1, WorkPerN
  if (WorkBasePerMeans(XReg,XPer).NE.MissVal) then		! calc anom v base
    PerTot = 0.0
    PerEn  = 0.0
    
!    write (99,"(2i6)"), XReg, XPer
    
    do XYear = 1, WorkYearN
      if (PertSelect(XYear).EQ.XPer) then
        if (WorkGotFull(XReg,XYear).NE.MissVal) then
          PerTot = PerTot + WorkGotFull(XReg,XYear)
          PerEn  = PerEn  + 1.0
        end if
      end if
    end do
    
!    write (99,"(2f10.4)"), PerTot, PerEn
    
    if (PerEn.GE.PertThreshold) then
      TotalValid = TotalValid + 1.0
      
      if (QAnomType.EQ.1) then
    	WorkAnomPerMeans(XReg,XPer) = (PerTot / PerEn) - WorkBasePerMeans(XReg,XPer)
      else if (QAnomType.EQ.2) then
    	WorkAnomPerMeans(XReg,XPer) = 100.0 * ((PerTot / PerEn) - WorkBasePerMeans(XReg,XPer)) / &
    					WorkBasePerMeans(XReg,XPer)
      end if
    end if

!    write (99,"(2f10.4)"), WorkBasePerMeans(XReg,XPer), WorkAnomPerMeans(XReg,XPer)
    
  end if
 end do
end do

PerCentValid = 100.0 * TotalValid / real(WorkRegN*WorkPerN)
print*, "  > Percentage period-region anoms valid: ", PerCentValid
    
deallocate (WorkBasePerMeans)

end subroutine CalcAnomMeans

!*******************************************************************************
! calc mean anom for globe for each period

subroutine GetWeights

print*, "  > Load region weights from .glo file ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), QGetWeights
	if (ReadStatus.LE.0.AND.QGetWeights.GE.1.AND.QGetWeights.LE.2) exit
end do

if (QGetWeights.EQ.1) then
  allocate ( Weights (WorkRegN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcPerGlo: Allocation failure #####"
  Weights = 1.0
else
  call LoadGlo (WorkLongN, WorkLatN, WorkRegN, WorkMapIDLReg, Weights)
end if

end subroutine GetWeights

!*******************************************************************************
! calc mean anom for globe for each period
! the mean anom is weighted by Weights()

subroutine CalcPerGlo

allocate ( PerGlo (WorkPerN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcPerGlo: Allocation failure #####"
PerGlo = MissVal

RegThresh = WorkRegN*(100.0-MissAccept)/100.0
PerValid  = 0.0

do XPer = 1, WorkPerN
  PerTot   = 0.0
  PerEn    = 0.0
  
  do XReg = 1, WorkRegN
    if (WorkAnomPerMeans(XReg,XPer).NE.MissVal) then
      PerTot = PerTot + (WorkAnomPerMeans(XReg,XPer) * Weights(XReg))
      PerEn  = PerEn  + Weights(XReg)
    end if
  end do

  if (PerEn.GE.RegThresh) then
  	PerGlo(XPer) = PerTot / PerEn
  	PerValid     = PerValid + 1.0
  end if
end do

print*, "  > Periods with valid weighted global means: ", PerValid

end subroutine CalcPerGlo

!*******************************************************************************
! calc JFB pattern noise (global)
! the weighting follows Mitchell et al 99 as far as possible
! the difference between local and global is squared BEFORE being weighted

subroutine CalcJFBGlo

NoiseJFB    = MissVal
NoiseJFBPop = MissVal

PerSqTot = 0.0
PerTot   = 0.0
PerEn    = 0.0

PerThresh = WorkPerN*(100.0-MissAccept)/100.0
RegThresh = WorkRegN*(100.0-MissAccept)/100.0

do XPer = 1, WorkPerN
 if (PerGlo(XPer).NE.MissVal) then
  RegSqTot    = 0.0
  RegEn       = 0.0

  do XReg = 1, WorkRegN
      if (WorkAnomPerMeans(XReg,XPer).NE.MissVal) then
        RegSqTot = RegSqTot + ( ((WorkAnomPerMeans(XReg,XPer) - PerGlo(XPer)) ** 2) * Weights(XReg) )
        RegEn    = RegEn    +    Weights(XReg)
      end if
  end do
  
  if (RegEn.GE.RegThresh) then
    RegSignal = RegSqTot / RegEn
    RegSignal = sqrt (RegSignal)
    write (99,"(i4,e16.5)"), XPer, RegSignal		! ###########
    
    PerSqTot  = PerSqTot + (RegSignal ** 2)
    PerTot    = PerTot   +  RegSignal
    PerEn     = PerEn    +  1.0
  end if
 end if
end do

if (PerEn.GE.PerThresh.AND.PerEn.GE.2) then
  NoiseJFB = PerTot / PerEn
  print*, "  > Mean  of signals from control periods: ", NoiseJFB

  NoiseJFB = (PerSqTot / PerEn) - ( (PerTot / PerEn) ** 2)
  NoiseJFB = (NoiseJFB * PerEn) / ( PerEn - 1 )
  NoiseJFB = sqrt (NoiseJFB) 
  print*, "  > Stdev of signals from control periods: ", NoiseJFB
else
  print*, "  > Too many missing signals to calc noise."
  print*, "  > Valid total and threshold values: ", PerEn, PerThresh
end if

end subroutine CalcJFBGlo

!*******************************************************************************
! calc stdev

subroutine CalcStDev

print*, "  > Calculate sample (=1) or population (=2) noise ?"
do
	read (*,*,iostat=ReadStatus), QSampPop
	if (ReadStatus.LE.0.AND.QSampPop.GE.1.AND.QSampPop.LE.2) exit
end do

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

TotalValid = 0.0
RMSEThreshold = WorkPerN*(100.0-MissAccept)/100.0
RegThresh     = WorkRegN*(100.0-MissAccept)/100.0
TotRMSE = 0.0
EnRMSE  = 0.0

write (99,"(A17,F10.4)"), "stdev threshold: ", RMSEThreshold

do XReg = 1, WorkRegN
  PerSqTot = 0.0
  PerTot   = 0.0
  PerEn    = 0.0
  
  do XPer = 1, WorkPerN
    if (WorkAnomPerMeans(XReg,XPer).NE.MissVal) then
      PerSqTot = PerSqTot + (WorkAnomPerMeans(XReg,XPer) ** 2)
      PerTot   = PerTot   +  WorkAnomPerMeans(XReg,XPer)
      PerEn    = PerEn    + 1
    end if
  end do
  
  if (PerEn.GE.RMSEThreshold) then
    if (QSampPop.EQ.2) then
      if (PerEn.GT.1) then
        AllRMSE (XReg) = (PerSqTot / PerEn) + ((PerTot / PerEn) * (PerTot / PerEn))
        AllRMSE (XReg) = AllRMSE (XReg) * PerEn / (PerEn - 1)
        AllRMSE (XReg) = sqrt (AllRMSE (XReg))
        TotalValid     = TotalValid + 1.0
      end if
    else
        AllRMSE (XReg) = (PerSqTot / PerEn) + ((PerTot / PerEn) * (PerTot / PerEn))
        AllRMSE (XReg) = sqrt (AllRMSE (XReg))
        TotalValid     = TotalValid + 1.0
    end if
    
    if (AllRMSE(XReg).GT.99999.99) AllRMSE(XReg) = 99999.99	! the .glo file cannot store larger values
    TotRMSE = TotRMSE + AllRMSE (XReg)
    EnRMSE  = EnRMSE  + 1.0    
  end if

!  write (99,"(i6,3f10.4)"), XReg, PerSqTot, PerTot, PerEn
end do

PerCentValid = 100.0 * TotalValid / real(WorkRegN)

print*, "  > Percentage region stdev valid:         ", PerCentValid

if (EnRMSE.GE.1.AND.EnRMSE.GT.RegThresh) &
	print*, "  > Global mean RMSE:                     ", TotRMSE/EnRMSE

print*, "  > Save the stdev to .glo file."
call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,AllRMSE,WorkMapIDLReg)

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

end subroutine CalcStDev

!*******************************************************************************
! calc stdev of stdev

subroutine CalcStDevStDev

print*, "  > How many members in ensemble to mimic ?"
do
	read (*,*,iostat=ReadStatus), EnsMemN
	if (ReadStatus.LE.0.AND.EnsMemN.GE.1.AND.EnsMemN.LT.WorkPerN) exit
end do

print*, "  > Calculate sample (=1) or population (=2) noise ?"
do
	read (*,*,iostat=ReadStatus), QSampPop
	if (ReadStatus.LE.0.AND.QSampPop.GE.1.AND.QSampPop.LE.2) exit
end do

print*, "  > Enter random number seed integer."
do
	read (*,*,iostat=ReadStatus), Seed (1)		! seed = rank-one size-one integer array
	if (ReadStatus.LE.0) exit
end do

WorkTrialN  = 200
MemThresh   = (100.0 - MissAccept) * real(EnsMemN) / 100.0 
TrialThresh = (100.0 - MissAccept) * WorkTrialN / 100.0 

allocate ( OpProc  (WorkTrialN,WorkRegN), &
           OpOut   (WorkTrialN),          &
           UsePer  (EnsMemN),             &
           LowerBound (WorkRegN),	  &
           UpperBound (WorkRegN),	  stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcStDevStDev: Allocation failure #####"
OpProc     = MissVal
OpOut      = MissVal
LowerBound = MissVal
UpperBound = MissVal

call random_seed (put=Seed)

do XTrial = 1, WorkTrialN
  UsePer = 0
  do XMem = 1, EnsMemN
    do
      call random_number (RandomReal)
      RandomReal    = (RandomReal * WorkPerN) + 0.5
      RandomMem     = nint (RandomReal)
      if (RandomMem.GE.1.AND.RandomMem.LE.WorkPerN) exit
    end do
    
    UsePer (XMem) = RandomMem
  end do
  
  write (99,*), XTrial, (UsePer(XMem),XMem=1,EnsMemN)	! #############
  
  do XReg = 1, WorkRegN
    PerSqTot = 0.0
    PerTot   = 0.0
    PerEn    = 0.0
    
    do XMem = 1, EnsMemN
      if (WorkAnomPerMeans(XReg,UsePer(XMem)).NE.MissVal) then
        PerSqTot = PerSqTot + (WorkAnomPerMeans(XReg,UsePer(XMem)) ** 2)
        PerTot   = PerTot   +  WorkAnomPerMeans(XReg,UsePer(XMem))
        PerEn    = PerEn    +  1.0
      end if
    end do
    
    if (PerEn.GE.MemThresh.AND.PerEn.GE.2) then
      OpProc (XTrial,XReg) = (PerSqTot / PerEn) - ( (PerTot / PerEn) ** 2 )
      
      if (QSampPop.EQ.2.AND.OpProc(XTrial,XReg).GT.0) &
      		OpProc (XTrial,XReg) = OpProc (XTrial,XReg) * PerEn / (PerEn - 1)
      if (OpProc(XTrial,XReg).GT.0) &
      		OpProc (XTrial,XReg) = sqrt (OpProc (XTrial,XReg))
    end if
  end do
end do

print*, "  > Is this a ref distribution for a one (=1) or two (=2) tailed test ?"
do
	read (*,*,iostat=ReadStatus), TailN
	if (ReadStatus.LE.0 .AND. TailN.GE.1 .AND. TailN.LE.2) exit
end do

print*, "  > What is the significance level (e.g. 0.95) ?"
do
	read (*,*,iostat=ReadStatus), SigLevel
	if (ReadStatus.LE.0 .AND. SigLevel.GE.0.5 .AND. SigLevel.LT.1) exit
end do

do XReg = 1, WorkRegN
  RegMiss = 0
  
  do XTrial = 1, WorkTrialN
    OpOut (XTrial) = OpProc (XReg,XTrial)
    if (OpProc (XReg,XTrial).EQ.MissVal) RegMiss = RegMiss + 1
  end do
  
  if ((WorkTrialN-RegMiss).GE.TrialThresh) then
    call BubbleSort (WorkTrialN,TotMiss,OpOut)

    CalcReal = (WorkTrialN - TotMiss) * (1.0 - SigLevel) / TailN
    CalcInt  = floor (CalcReal)
    LowerBound (XReg) = OpOut (CalcInt)
  
    CalcReal = (WorkTrialN - TotMiss) * SigLevel / TailN
    CalcInt  = ceiling (CalcReal)
    UpperBound (XReg) = OpOut (CalcInt)  
  end if

  write (99,"(i6,2f10.4)"), XReg, LowerBound (XReg), UpperBound (XReg)
end do

if (TailN.EQ.1) print*, "  > Save sig value for 1-tail (lower): "
if (TailN.EQ.2) print*, "  > Save sig value for 2-tail (lower): "
call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,LowerBound,WorkMapIDLReg)

if (TailN.EQ.1) print*, "  > Save sig value for 1-tail (upper): "
if (TailN.EQ.2) print*, "  > Save sig value for 2-tail (upper): "
call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,UpperBound,WorkMapIDLReg)

deallocate (LowerBound,UpperBound,OpOut,UsePer,OpProc, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcStDevStDev: Deallocation failure #####"

end subroutine CalcStDevStDev

!*******************************************************************************
! calc interann var for ens of 30y slices: save mean+sd of set
! no anomalising, so full set of PerLen+GapLen within control
! that full set is divided by EnsMemN, to give the sample size
! if PerLen=30, GapLen=10, EnsMemN=4, and WorkYearN=1400, then
!    sample of 8:
!  i = periods 1, 9,17,25
! ii = periods 2,10,18,26 etc
! from the set of seasonal means (120 per sample member in e.g.), we calc interann var
! for the sample, we calc mean and stdev-pop

subroutine InterAnnVar

print*, "  > The variance is the same whether or not we anom(abs), so we don't."
call Specification1		! gives PerLen and GapLen
call IdentifyPer3		! gives WorkPerN and PertSelect(XYear)

print*, "  > How many members in ensemble to mimic ?"
do
	read (*,*,iostat=ReadStatus), EnsMemN
	if (ReadStatus.LE.0.AND.EnsMemN.GE.1.AND.EnsMemN.LT.WorkPerN) exit
end do
EnsPerYearN = EnsMemN * PerLen
EnsPerN     = floor (real(WorkPerN)/real(EnsMemN))

allocate ( OpOut  (EnsPerYearN), 	&
	   OpProc (WorkRegN,EnsPerN), 	stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: InterAnnVar: Allocation failure #####"
OpProc = MissVal

do XReg = 1, WorkRegN
  do XPer = 1, EnsPerN
    PerYear = 0
    OpOut   = MissVal
    
    do XMem = 1, EnsMemN
     do XYear = 1, WorkYearN    

      ThisPer = XPer + ( (XMem-1) * EnsMemN ) 				! select per to include

      if (PertSelect(XYear).EQ.ThisPer) then
      	PerYear = PerYear + 1
      	OpOut(PerYear) = WorkGotFull (XReg,XYear)			! fill vector with PerLen seasonals
      end if

     end do
    end do
    
    OpProc (XReg,XPer) = OpCalcStDev (OpOut, EnsPerYearN, MissAccept, 3) ! calc samp-var (=3) for Reg and Per 
  end do
end do

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

allocate ( OpOut  (WorkPerN), &
	   AllAye (WorkRegN), &
	   AllBee (WorkRegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: InterAnnVar: Allocation failure #####"

do XReg = 1, WorkRegN
  do XPer = 1, EnsPerN
    OpOut (XPer) = OpProc (XReg,XPer)					! fill vector with variances
  end do
  
  AllAye (XReg) = OpCalcMean  (OpOut, EnsPerN, MissAccept)		! calc samp-mean for Reg
  AllBee (XReg) = OpCalcStDev (OpOut, EnsPerN, MissAccept, 1)		! calc pop-sd (=2) for Reg
end do

print*, "  > Save means of period variances to .glo:"
call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,AllAye,WorkMapIDLReg)

print*, "  > Save sd-pop of period variances to .glo:"
call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,AllBee,WorkMapIDLReg)

deallocate (OpOut,OpProc,AllAye,AllBee, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: InterAnnVar: Deallocation failure #####"

end subroutine InterAnnVar

!*******************************************************************************
! end

end program SamplingError2