! ScalePattern.f90
! written by Tim Mitchell on 07.11.01
! last modified on 07.11.01
! program to scale response pattern from grip file
! designed using the linked list concept
! f90 -o ./../grim/scalepattern sortmod.f90 filenames.f90 initialmod.f90 time.f90 grimfiles.f90 annfiles.f90 
!	glofiles.f90 aggfiles.f90 perfiles.f90 regress.f90 ./../grim/scalepattern.f90

program ScalePattern

use FileNames
use InitialMod
use Time
use GrimFiles
use AnnFiles
use GloFiles
use AggFiles
use PerFiles
use Regress

implicit none

type Exec
  character (len=20) 				:: Name			! execution name
  type (Exec), pointer 				:: Prev			! prev execution defined: recursion
  type (Exec), pointer 				:: Next			! next execution defined: recursion

  character (len=80)				:: Kay1File, Kay2File	! pattern grip files
  character (len=80)				:: GloTFile, EquTFile	! .ann for target
  character (len=80)				:: TargFile, BaseFile	! target grim files
end type Exec

type (Exec), pointer				:: OneExec, CurrentExec, StackExec

character (len=80), parameter 	:: SpecFile = "/cru/mikeh1/f709762/f90/grim/_ref/ops/scalepattern.21.ops.X"

real, pointer, dimension (:,:,:)		:: Est,Mod
real, pointer, dimension (:,:,:)		:: Array3D
real, pointer, dimension (:,:)			:: Array2D,ArrayUM,ArrayWM,ArrayUA,ArrayWA
real, pointer, dimension (:,:)			:: Kay1,Kay2 		! the arrays of k: SeasonN, BoxN
real, pointer, dimension (:,:)			:: GloWeights
real, pointer, dimension (:)			:: GloT,EquT,CurrentTS  ! ts: global T, equilibrium T
real, pointer, dimension (:)			:: Array1D,TSLowVec,TSHighVec,Extract
real, pointer, dimension (:)			:: Weight1D,RegWeights

integer, pointer, dimension (:,:)		:: Grid, LoadGrid, MapIDLReg
integer, pointer, dimension (:)			:: YearAD, LoadYearAD, RegSizes
integer, pointer, dimension (:)			:: DumpYear		! periods to dump

character (len=80), pointer, dimension (:)	:: DumpRegSet,TextRegSet,RegNamesUniq80,SaveRegFile
character (len=20), pointer, dimension (:) 	:: RegNames,RegNamesUniq20
character (len= 9),  pointer, dimension(:)	:: ColTitles

real, dimension (4)				:: Bounds, LoadBounds

integer, dimension (3)				:: DumpStat
integer, dimension (17)				:: DumpCal

character (len=3), dimension (17), parameter	:: SeasonNames = ['jan','feb','mar','apr','may','jun','jul',&
						'aug','sep','oct','nov','dec','MAM','JJA','SON','DJF','ann']

character (len=3), dimension (3), parameter	:: StatsText = ['est','mod','dif']

real, parameter 	:: MissVal = -999.0
integer, parameter 	:: SeasonN=17, MonthN=12, BoundN=4, StatN=3

character (len=80), parameter 	:: Blank    = ""
real :: MissAccept, MissThresh, OpTot, OpEn, OpMiss, VariFactor
real :: OpTotUM,OpTotUA,OpTotWM,OpTotWA,OpEnWei,OpMissWei

integer :: ReadStatus, AllocStat
integer :: ExecN, BoxN, ExeN, WyeN, YearN, LoadYearN, DumpYearN, DumpRegSetN, RegN
integer :: XExec, XBox, XExe, XWye, XYear, XLoadYear, XDumpYear, XDumpRegSet, XReg
integer :: XSeason, XMonth, XBound, XStat
integer :: PerLen, GapLen, YearAD0, YearAD1, YearLimit, CheckBoxN, ThisYear, ThisMonth
integer :: QMethod,QCompare,QMeanSum,QDumpGrip,QDumpGlo,QDumpAgg,QDumpGrim,QDumpPer,QDumpGlobe
integer :: QGloAgg,QPerGlobe
integer :: Year0,Year1,LoadYear0,LoadYear1,CallVariCode

character (len=80) :: LoadInfo,LoadFile,SaveInfo,SaveFile,GridRef,VariName
character (len=80) :: TextDir,TextInfo,TextTarg,TextPatt,SaveSTem,SaveTip
character (len=20) :: TextYear
character (len= 4) :: LoadSuffix,SaveSuffix,Variable

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

call Intro
call GetOpsFile
call SetSpecifics

call ResetCurrentExec		! points CurrentExec back to the first exec

do
  print*, "  > Execution: ", CurrentExec%Name
  
  Variable = trim(CurrentExec%Name)
  if (Variable.EQ.".wet".OR.Variable.EQ.".pre".OR.Variable.EQ.".frs") then
    QMeanSum = 2
  else
    QMeanSum = 1
  end if
  
  print*, "  > Loading patterns..."
  call LoadPatterns  
  print*, "  > Loading scalers..."
  call LoadScalers
  print*, "  > Making estimate..."
  call MakeEstimate
  
  if (QCompare.EQ.1) then
    print*, "  > Making comparison..."
    call LoadModelled
!    call CompareEstMod			! has to be omitted because of lack of memory - workarounds used
  end if 
  
  print*, "  > Dumping to file..."
  if (QDumpGrip .EQ.1) call DumpGrip  
  if (QDumpGlo  .EQ.1) then
  	QGloAgg = 1 ; call DumpGloAgg
  end if
  if (QDumpAgg  .EQ.1) then
  	QGloAgg = 2 ; call DumpGloAgg
  end if
  if (QDumpGrim .EQ.1) call DumpGrim
  if (QDumpPer  .EQ.1) call DumpPer
  if (QDumpGlobe.EQ.1) call DumpGlobe
  
  call RemoveEstMod
  
  if (.not.associated(CurrentExec%Next)) exit					! last exec, so exit loop
  if (associated(CurrentExec%Next)) CurrentExec => CurrentExec%Next		! not 1st exec, so next exec  
end do

call Finish

contains

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

subroutine Intro

open (99,file="/cru/mikeh1/f709762/scratch/log/log-scalepattern.dat",status="replace",action="write")
print*
print*, "  > ***** ScalePattern: scales response patterns *****"
print*
print*, "  > Spec file: ", trim(SpecFile)
print*

end subroutine Intro

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

subroutine GetOpsFile

open  (1,file=SpecFile,status="old",access="sequential",form="unformatted",action="read")

read  (1), ExecN, BoxN, ExeN, WyeN
read  (1), PerLen, MissAccept, QMethod, QCompare, YearAD0, YearAD1
read  (1), (Bounds(XBound),XBound=1,4)
read  (1), (DumpStat(XStat),XStat=1,3)
read  (1), (DumpCal(XSeason),XSeason=1,17)
read  (1), DumpYearN,DumpRegSetN
read  (1), QDumpGrip,QDumpGlo,QDumpAgg,QDumpGrim,QDumpPer,QDumpGlobe
read  (1), TextInfo,TextDir,TextTarg,TextPatt

if (DumpYearN.GT.0) then 
  allocate (DumpYear (DumpYearN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Intro: Allocation failure: DumpYear #####"

  read  (1), (DumpYear(XDumpYear),XDumpYear=1,DumpYearN)
end if

if (DumpRegSetN.GT.0) then 
  allocate (DumpRegSet (DumpRegSetN), &
  	    TextRegSet (DumpRegSetN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Intro: Allocation failure: DumpRegSet #####"

  read  (1), (DumpRegSet(XDumpRegSet),XDumpRegSet=1,DumpRegSetN)
  read  (1), (TextRegSet(XDumpRegSet),XDumpRegSet=1,DumpRegSetN)
end if

allocate (Grid   (ExeN,WyeN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Intro: Allocation failure: Grid #####"
do XExe = 1, ExeN
  read (1), (Grid(XExe,XWye), XWye=1,WyeN)
end do

do XExec = 1, ExecN
  allocate (OneExec, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Intro: Allocation failure: OneExec #####"
  
  read  (1), OneExec%Name
  read  (1), OneExec%Kay1File, OneExec%Kay2File, OneExec%GloTFile, OneExec%EquTFile
  read  (1), OneExec%BaseFile, OneExec%TargFile
  
  if      (XExec.EQ.1) then
    nullify(OneExec%Prev)
    nullify(OneExec%Next)
  else
    OneExec%Prev     => CurrentExec
    CurrentExec%Next => OneExec
    nullify(OneExec%Next)
  end if
  
  CurrentExec => OneExec
end do
  
close (1)

end subroutine GetOpsFile

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

subroutine SetSpecifics

YearN = YearAD1 - YearAD0 + 1

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

do XYear = 1, YearN
    YearAD (XYear) = YearAD0 + XYear - 1
end do

GridRef = GetGloRef (ExeN,WyeN)

end subroutine SetSpecifics

!*******************************************************************************
! kay1 (and kay2 if required) is in a grip file
! it is expressed in <variable units> per degC

subroutine LoadPatterns

if (CurrentExec%Kay1File.NE."") then
  allocate (Kay1(SeasonN,BoxN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadPatterns: Allocation failure: Kay1 #####"
  Kay1 = MissVal
  
  LoadFile = CurrentExec%Kay1File
  Array2D => Kay1
  call GrabGrip
end if
 
if (CurrentExec%Kay2File.NE."") then
  allocate (Kay2(SeasonN,BoxN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadPatterns: Allocation failure: Kay2 #####"
  Kay2 = MissVal

  LoadFile = CurrentExec%Kay2File
  Array2D => Kay2
  call GrabGrip
end if

nullify (Array2D)

end subroutine LoadPatterns

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

subroutine GrabGrip

call LoadGrip (Array2D,LoadGrid,LoadBounds,LoadInfo,LoadFile,CurrentExec%Name,LoadSuffix)

call CheckGridAB (Grid,LoadGrid,CheckBoxN)
if (CheckBoxN.EQ.MissVal) then
  print*, "  > ##### ERROR: GrabGrip: grids do not match #####"
  
  deallocate (Array2D,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabGrip: Deallocation failure #####" 
end if

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

end subroutine GrabGrip

!*******************************************************************************
! GloT (and EquT if required) should be in a single-column .ann file
! both should already be anomalised relative to the same base period as in the BaseFile
! both should already be smoothed at the periodicity PerLen

subroutine LoadScalers

if (CurrentExec%GloTFile.NE."") then
  allocate (GloT(YearN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadScalers: Allocation failure: GloT #####"
  GloT = MissVal

  call LoadAnn (CurrentExec%GloTFile, LoadYearAD, ColTitles, Array2D)
  CurrentTS => GloT
  call GrabAnn
end if

if (CurrentExec%EquTFile.NE."") then
  allocate (EquT(YearN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadScalers: Allocation failure: EquT #####"
  EquT = MissVal

  call LoadAnn (CurrentExec%EquTFile, LoadYearAD, ColTitles, Array2D)
  CurrentTS => EquT
  call GrabAnn
end if

nullify (CurrentTS)

end subroutine LoadScalers

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

subroutine GrabAnn

call CommonVecPer (YearAD,LoadYearAD,Year0,Year1,LoadYear0,LoadYear1)

if (size(Array2D,2).EQ.1) then
 if (Year0.NE.MissVal) then
  XYear = Year0 - 1 ; XLoadYear = LoadYear0 - 1
  do
    XYear = XYear + 1 ; XLoadYear = XLoadYear + 1    
    CurrentTS (XYear) = Array2D (XLoadYear,1)
    if (XYear.EQ.Year1) exit
  end do
 else
  print*, "  > ##### ERROR: GrabAnn: no period in common #####"
 end if
else
 print*, "  > ##### ERROR: GrabAnn: multiple columns in .ann #####"
end if

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

end subroutine GrabAnn

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

subroutine MakeEstimate

allocate (Est(YearN,SeasonN,BoxN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: MakeEstimate: Allocation failure: Est #####"
Est = MissVal

if      (QMethod.EQ. 1) then
  do XYear = 1, YearN
    do XSeason = 1, SeasonN
      do XBox = 1, BoxN
        Est(XYear,XSeason,XBox) = Kay1(XSeason,XBox) * GloT(XYear)
      end do
    end do
  end do

  deallocate (Kay1,GloT, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: MakeEstimate: Deallocation failure: 01 #####"
else if (QMethod.EQ.11) then
  do XYear = 1, YearN
    do XSeason = 1, SeasonN
      do XBox = 1, BoxN
        Est(XYear,XSeason,XBox) = (Kay1(XSeason,XBox) * GloT(XYear)) + &
        				(Kay2(XSeason,XBox) * (EquT(XYear) - GloT(XYear)))
      end do
    end do
  end do

  deallocate (Kay1,Kay2,GloT,EquT, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: MakeEstimate: Deallocation failure: 11 #####"
else
  print*, "  > ##### ERROR: MakeEstimate: QMethod unrecognised #####"
end if

end subroutine MakeEstimate

!*******************************************************************************
! we take the unanomalised, unsmoothed, target that we are aiming at, as a grim file
! we anomalise the target versus the base file (which must be a single year)
! we fill in the seasonal (cols 13-16) and annual (17) means
! we smooth the time series at each box and period (1...17) at periodicity PerLen

subroutine LoadModelled

allocate (Mod(YearN,SeasonN,BoxN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadModelled: Allocation failure: Mod #####"
Mod = MissVal

call LoadTarg
call LoadBase
call GetMod1317
call SmooMod

end subroutine LoadModelled

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

subroutine LoadTarg

call LoadGrim (Array3D,LoadGrid,LoadYearAD,LoadBounds,LoadInfo,CurrentExec%TargFile,&
		CurrentExec%Name,LoadSuffix)

call CommonVecPer (YearAD,LoadYearAD,Year0,Year1,LoadYear0,LoadYear1)

call CheckGridAB (Grid,LoadGrid,CheckBoxN)

if (Year0.EQ.MissVal) then
 print*, "  > ##### ERROR: LoadTarg: no period in target file that matches YearAD #####"
else
 if (CheckBoxN.EQ.MissVal) then
  print*, "  > ##### ERROR: LoadTarg: grids do not match #####"
 else
  do XMonth = 1, MonthN
   do XBox = 1, BoxN
    do XYear = Year0, Year1
     XLoadYear = LoadYear0 + XYear - Year0     
     if (Array3D(XLoadYear,XMonth,XBox).NE.MissVal) Mod (XYear,XMonth,XBox) = Array3D (XLoadYear,XMonth,XBox)
    end do
   end do
  end do
 end if
end if

deallocate (LoadGrid,LoadYearAD,Array3D,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadTarg: : Deallocation failure #####"

end subroutine LoadTarg

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

subroutine LoadBase

call LoadGrim (Array3D,LoadGrid,LoadYearAD,LoadBounds,LoadInfo,CurrentExec%BaseFile,&
		CurrentExec%Name,LoadSuffix)

call CheckGridAB (Grid,LoadGrid,CheckBoxN)

if (size(LoadYearAD).NE.1) then
 print*, "  > ##### ERROR: LoadBase: grim not a single year #####"
else
 if (CheckBoxN.EQ.MissVal) then
  print*, "  > ##### ERROR: LoadBase: grids do not match #####"
 else
  do XMonth = 1, MonthN
   do XBox = 1, BoxN
    if (Array3D(1,XMonth,XBox).NE.MissVal) then
      do XYear = 1, YearN
        if (Mod(XYear,XMonth,XBox).NE.MissVal) Mod (XYear,XMonth,XBox) = &
        				Mod (XYear,XMonth,XBox) - Array3D (1,XMonth,XBox)
      end do
    else
      do XYear = 1, YearN
        Mod (XYear,XMonth,XBox) = MissVal
      end do
    end if
   end do
  end do
 end if
end if

deallocate (LoadGrid,LoadYearAD,Array3D,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadBase: : Deallocation failure #####"

end subroutine LoadBase

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

subroutine GetMod1317

do XYear = 1, YearN				! get annual data
  do XBox = 1, BoxN
    OpTot = 0.0 ; OpEn = 0.0
    
    do XMonth = 1, MonthN
      if (Mod(XYear,XMonth,XBox).NE.MissVal) then
        OpTot = OpTot + Mod(XYear,XMonth,XBox)
        OpEn  = OpEn  + 1.0
      end if
    end do
    
    if (OpEn.EQ.12) then
    	if (QMeanSum.EQ.1) Mod(XYear,17,XBox) = OpTot / OpEn
    	if (QMeanSum.EQ.2) Mod(XYear,17,XBox) = OpTot
    end if
  end do
end do

do XYear = 1, YearN				! get seasonal data
 do XBox = 1, BoxN
  do XSeason = 13, 16       
   OpTot = 0.0 ; OpEn = 0.0
    
   do XMonth = (((XSeason-13)*3)+3), (((XSeason-13)*3)+5)
      ThisYear = XYear ; ThisMonth = XMonth
      if (XMonth.GT.12) then
        ThisYear = XYear + 1 ; ThisMonth = XMonth - 12
      end if
      
      if (ThisYear.LE.YearN) then
       if (Mod(ThisYear,ThisMonth,XBox).NE.MissVal) then
        OpTot = OpTot + Mod(ThisYear,ThisMonth,XBox)
        OpEn  = OpEn  + 1.0
       end if
      end if
   end do
    
   if (OpEn.EQ.3) then
    	if (QMeanSum.EQ.1) Mod(XYear,XSeason,XBox) = OpTot / OpEn
    	if (QMeanSum.EQ.2) Mod(XYear,XSeason,XBox) = OpTot
   end if
  end do
 end do
end do

end subroutine GetMod1317

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

subroutine SmooMod

allocate (Array1D  (YearN), &
	  TSLowVec (YearN), &
	  TSHighVec(YearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SmooMod: Allocation failure #####"

do XSeason = 1, SeasonN
  do XBox = 1, BoxN
    TSLowVec = MissVal ; TSHighVec = MissVal
    
    do XYear = 1, YearN
      Array1D (XYear) = Mod (XYear,XSeason,XBox)
    end do
    
    call GaussSmooth (YearN,PerLen,1,Array1D,TSLowVec,TSHighVec)
    
    do XYear = 1, YearN
      Mod (XYear,XSeason,XBox) = TSLowVec (XYear)
    end do
  end do
end do
    
deallocate (Array1D,TSLowVec,TSHighVec, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SmooMod: Deallocation failure: Array1D #####"

end subroutine SmooMod

!*******************************************************************************
! Dif = Est -- Mod
! has to be omitted and worked around because of lack of memory
!
!subroutine CompareEstMod
!
!allocate (Dif(YearN,SeasonN,BoxN), stat=AllocStat)
!if (AllocStat.NE.0) print*, "  > ##### ERROR: CompareEstMod: Allocation failure: Dif #####"
!Dif = MissVal
!
!do XYear = 1, YearN
!  do XSeason = 1, SeasonN
!    do XBox = 1, BoxN
!      if (Mod(XYear,XSeason,XBox).NE.MissVal.AND.Est(XYear,XSeason,XBox).NE.MissVal) then
!        Dif(XYear,XSeason,XBox) = Est(XYear,XSeason,XBox) - Mod(XYear,XSeason,XBox)
!      end if
!    end do
!  end do
!end do
!
!end subroutine CompareEstMod
!
!*******************************************************************************

subroutine DumpGrip

write (99,*), "Dumping to grip..."

allocate (Array2D(SeasonN,BoxN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpGrip: Allocation failure: Array2D #####"

do XDumpYear = 1, DumpYearN
  TextYear = GetTextFromInt (YearAD(DumpYear(XDumpYear)))
  
  do XStat = 1, StatN   
    if (DumpStat(XStat).EQ.1) then
      if      (XStat.EQ.1) then
       do XSeason = 1, SeasonN
        do XBox = 1, BoxN
          Array2D (XSeason,XBox) = Est (DumpYear(XDumpYear),XSeason,XBox)
        end do
       end do
      else if (XStat.EQ.2) then
       do XSeason = 1, SeasonN
        do XBox = 1, BoxN
          Array2D (XSeason,XBox) = Mod (DumpYear(XDumpYear),XSeason,XBox)
        end do
       end do
      else if (XStat.EQ.3) then
       do XSeason = 1, SeasonN
        do XBox = 1, BoxN
         if (Est(DumpYear(XDumpYear),XSeason,XBox).NE.MissVal.AND.Mod(DumpYear(XDumpYear),XSeason,XBox) &
         		.NE.MissVal) then
          Array2D (XSeason,XBox) = Est(DumpYear(XDumpYear),XSeason,XBox)-Mod(DumpYear(XDumpYear),XSeason,XBox)
         else
          Array2D (XSeason,XBox) = MissVal
         end if
        end do
       end do
      end if
      
      SaveInfo = trim(TextInfo) // " " // trim(TextTarg) // "@" // trim(TextYear) // &
      			" stat=" // StatsText(XStat)       				
      if (XStat.NE.2) SaveInfo = trim(SaveInfo) // " patt=" // trim(TextPatt)
      if (QMeanSum.EQ.1) SaveInfo = trim(SaveInfo) // " time=mean"
      if (QMeanSum.EQ.2) SaveInfo = trim(SaveInfo) // " time=sum"
      
      SaveFile = trim(TextDir) // StatsText(XStat) // "-" // trim(TextTarg) // "-" // trim(TextYear)
      if (XStat.NE.2) SaveFile = trim(SaveFile) // "." // trim(TextPatt)
      SaveFile = trim(SaveFile) // CurrentExec%Name
   
      call SaveGrip (Array2D,Grid,Bounds,SaveInfo,SaveFile,CurrentExec%Name,SaveSuffix)
    end if
  end do
end do

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

end subroutine DumpGrip

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

subroutine DumpGloAgg

if (QGloAgg.EQ.1) write (99,*), "Dumping to .glo..."
if (QGloAgg.EQ.2) write (99,*), "Dumping to .agg..."

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

do XDumpRegSet = 1, DumpRegSetN
 call PrepareRegions
 
 allocate (Array1D   (RegN), &
 	   Weight1D  (RegN), stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpGloAgg: Allocation failure: Array1D #####"
 
 if (QGloAgg.EQ.2) then
   allocate (Array2D        (RegN,SeasonN), &
   	     RegNamesUniq20 (RegN),         stat=AllocStat)
   if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpGloAgg: Allocation failure: Array2D #####"   
   
   call MakeBatch (Blank,Blank,RegNames,RegNamesUniq80)
   RegNamesUniq20 = RegNamesUniq80
 end if
 
 do XStat = 1, StatN   
  if (DumpStat(XStat).EQ.1) then
    do XDumpYear = 1, DumpYearN
      if (QGloAgg.EQ.2) Array2D = 0.0
      
      do XSeason = 1, SeasonN
       if (DumpCal(XSeason).EQ.1.OR.QGloAgg.EQ.2) then		! checks whether this dump is required
        XYear = DumpYear(XDumpYear)
        call MakeBoxVector      
        call MakeRegVector
        
        if      (QGloAgg.EQ.1) then
          TextYear = GetTextFromInt (YearAD(DumpYear(XDumpYear)))
  
          SaveInfo = trim(TextInfo) // " " // trim(TextTarg) // "@" // trim(TextYear) // " stat=" // &
      			StatsText(XStat)       				
          if (XStat.NE.2) SaveInfo = trim(SaveInfo) // " patt=" // trim(TextPatt)
          SaveInfo = trim(SaveInfo) // " wei-regs=" // trim(TextRegSet(XDumpRegSet)) // " " // & 
      			SeasonNames(XSeason) // "-"
          if (QMeanSum.EQ.1) SaveInfo = trim(SaveInfo) // "mean"
          if (QMeanSum.EQ.2) SaveInfo = trim(SaveInfo) // "sum"
      
          SaveFile = trim(TextDir) // StatsText(XStat) // "-" // trim(TextTarg) // "-" // trim(TextYear)
          if (XStat.NE.2) SaveFile = trim(SaveFile) // "." // trim(TextPatt)
          SaveFile = trim(SaveFile) // "." // trim(TextRegSet(XDumpRegSet)) // "." // SeasonNames(XSeason) //&
      			trim(CurrentExec%Name) // ".glo"

          call SaveGlo (ExeN,WyeN,RegN,GridRef,SaveFile,SaveInfo,Array1D,MapIDLReg)
        else if (QGloAgg.EQ.2) then
	  do XReg = 1, RegN
            Array2D (XReg,XSeason) = Array1D(XReg)
          end do
        end if        
      end if
     end do
     
     if (QGloAgg.EQ.2) then
       TextYear = GetTextFromInt (YearAD(DumpYear(XDumpYear)))
       
       SaveInfo = trim(TextInfo) // " " // trim(TextTarg) // "@" // trim(TextYear) // " stat=" // &
      			StatsText(XStat)       				
       if (XStat.NE.2) SaveInfo = trim(SaveInfo) // " patt=" // trim(TextPatt)
       SaveInfo = trim(SaveInfo) // " wei-regs=" // trim(TextRegSet(XDumpRegSet)) // " "
       if (QMeanSum.EQ.1) SaveInfo = trim(SaveInfo) // "mean"
       if (QMeanSum.EQ.2) SaveInfo = trim(SaveInfo) // "sum"
      
       SaveFile = trim(TextDir) // StatsText(XStat) // "-" // trim(TextTarg) // "-" // trim(TextYear)
       if (XStat.NE.2) SaveFile = trim(SaveFile) // "." // trim(TextPatt)
       SaveFile = trim(SaveFile) // "." // trim(TextRegSet(XDumpRegSet)) // trim(CurrentExec%Name) // ".agg"

       call SaveAgg(SaveFile,RegNamesUniq20,AllData=Array2D,NoResponse=1)
     end if
   end do
  end if
 end do

 if (QGloAgg.EQ.2) then
   deallocate (Array2D,RegNamesUniq20,RegNamesUniq80, stat=AllocStat)
   if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpGloAgg: Deallocation failure: Array2D #####"
 end if
 
 deallocate (Array1D,Weight1D,GloWeights,RegWeights,MapIDLReg,RegSizes,RegNames, stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpGloAgg: Deallocation failure: Array1D #####"
end do

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

end subroutine DumpGloAgg

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

subroutine DumpGrim

write (99,*), "Dumping ts at every box to grim..."

do XStat = 1, StatN   
  if (DumpStat(XStat).EQ.1) then
    SaveInfo = trim(TextInfo) // " " // trim(TextTarg) // "@" // trim(TextYear) // " stat=" // &
      			StatsText(XStat)       				
    if (XStat.NE.2) SaveInfo = trim(SaveInfo) // " patt=" // trim(TextPatt) // " "
    if (QMeanSum.EQ.1) SaveInfo = trim(SaveInfo) // "mean"
    if (QMeanSum.EQ.2) SaveInfo = trim(SaveInfo) // "sum"
      
    SaveFile = trim(TextDir) // StatsText(XStat) // "-" // trim(TextTarg)
    if (XStat.NE.2) SaveFile = trim(SaveFile) // "." // trim(TextPatt) // "." // trim(CurrentExec%Name)

    if (XStat.EQ.1) call SaveGrim (Est,Grid,YearAD,Bounds,SaveInfo,SaveFile,trim(CurrentExec%Name),SaveSuffix)
    if (XStat.EQ.2) call SaveGrim (Mod,Grid,YearAD,Bounds,SaveInfo,SaveFile,trim(CurrentExec%Name),SaveSuffix)

    if (XStat.EQ.3) then
      allocate (Array3D(YearN,SeasonN,BoxN), stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpGrim: Allocation failure: Array3D #####"      
      Array3D = MissVal
      
      do XYear = 1, YearN
        do XSeason = 1, SeasonN
          do XBox= 1, BoxN
            if (Est(XYear,XSeason,XBox).NE.MissVal.AND.Mod(XYear,XSeason,XBox).NE.MissVal) &
            		Array3D(XYear,XSeason,XBox) = Est(XYear,XSeason,XBox) - Mod(XYear,XSeason,XBox)
          end do
        end do
      end do
      
      call SaveGrim (Array3D,Grid,YearAD,Bounds,SaveInfo,SaveFile,trim(CurrentExec%Name),SaveSuffix)
      
      deallocate (Array3D, stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpGrim: Deallocation failure: Array3D #####"      
    end if
  end if
end do
  
end subroutine DumpGrim

!*******************************************************************************
! each seasonal and annual mean is calc as a mean (QMeanSum=1) or sum (=2)
! each regional value is calculated as the weighted mean of its constiuent boxes
! each region has its yearly/17 values dumped to .per
! this option should not be chosen except where the number of regions is relatively small

subroutine DumpPer

if (QMeanSum.EQ.1) CallVariCode = 0
if (QMeanSum.EQ.2) CallVariCode = 2
   
do XDumpRegSet = 1, DumpRegSetN
 call PrepareRegions
 
   allocate (Array2D        (RegN,SeasonN), &
   	     RegNamesUniq20 (RegN),         stat=AllocStat)
   if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpGloAgg: Allocation failure: Array2D #####"   
   
   call MakeBatch (Blank,Blank,RegNames,RegNamesUniq80)
   RegNamesUniq20 = RegNamesUniq80

 allocate (Array3D     (RegN,YearN,SeasonN), &
 	   Array2D     (YearN,SeasonN),      &
 	   Weight1D    (RegN),               &
 	   Array1D     (RegN),               &
 	   Extract     (BoxN),               &
 	   SaveRegFile (RegN),               stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpPer: Allocation failure: Array3D #####"
 
 SaveStem = trim(TextDir) // StatsText(XStat) // "-" // trim(TextTarg) // "."
 if (XStat.NE.2) SaveStem = trim(SaveStem) // trim(TextPatt) // "."
 SaveTip = "." // trim(CurrentExec%Name) // ".per"
  
 call MakeBatch (SaveStem,SaveTip,RegNames,SaveRegFile)
   
 do XStat = 1, StatN   
  if (DumpStat(XStat).EQ.1) then
   do XYear = 1, YearN
     do XSeason = 1, SeasonN
       call MakeBoxVector
       call MakeRegVector
       
       do XReg = 1, RegN
         Array3D (XReg,XYear,XSeason) = Array1D (XReg)
       end do 
     end do
   end do
   
   do XReg = 1, RegN
     do XYear = 1, YearN
       do XSeason = 1, SeasonN
         Array2D(XYear,XSeason) = Array3D(XReg,XYear,XSeason)
       end do
     end do
     
     call SavePer (SaveRegFile(XReg), YearAD, CallVariCode, AllData=Array2D, NoResponse=1)
   end do
  end if
 end do
 
 deallocate (MapIDLReg,RegSizes,RegNames,GloWeights,RegWeights,&
   		Array3D,Array2D,Array1D,Weight1D,Extract, stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpPer: Deallocation failure: Reg* #####"
end do

end subroutine DumpPer

!*******************************************************************************
! each seasonal and annual mean is calc as a mean (QMeanSum=1) or sum (=2)
! each regional value is calculated as the weighted mean of its constiuent boxes
! um: the global mean of all regions is unweighted
! ua: the global average absolute value of all regions is unweighted
! wm: the global mean of all regions is weighted
! wa: the global average absolute value of all regions is weighted

subroutine DumpGlobe

if (QMeanSum.EQ.1) CallVariCode = 0
if (QMeanSum.EQ.2) CallVariCode = 2
   
do XDumpRegSet = 1, DumpRegSetN
 call PrepareRegions
   
 allocate (ArrayUM   (YearN,SeasonN), &
           ArrayUA   (YearN,SeasonN), &
           ArrayWM   (YearN,SeasonN), &
           ArrayWA   (YearN,SeasonN), &
 	   Weight1D  (RegN),          &
 	   Array1D   (RegN),          &
 	   Extract   (BoxN),          stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpGlobe: Allocation failure: Array2D #####"
   
 do XStat = 1, StatN   
  if (DumpStat(XStat).EQ.1) then
   ArrayUM = MissVal ; ArrayUA = MissVal ; ArrayWM = MissVal ; ArrayWA = MissVal
   
   do XYear = 1, YearN
     do XSeason = 1, SeasonN
       call MakeBoxVector
       call MakeRegVector
       
       OpTotUM = 0.0 ; OpTotUA = 0 ; OpTotWM = 0.0 ; OpTotWA = 0
       OpEn = 0.0 ; OpEnWei = 0.0 ; OpMiss = 0.0 ; OpMissWei = 0.0
       do XReg = 1, RegN
         if (Array1D(XReg).NE.MissVal) then
           OpTotUM = OpTotUM + Array1D(XReg)
           OpTotUA = OpTotUA + abs(Array1D(XReg))
           OpEn    = OpEn    + 1.0
           
           if (Weight1D(XReg).NE.MissVal) then
             OpTotWM = OpTotWM + (Array1D(XReg)*Weight1D(XReg))
             OpTotWA = OpTotWA + (abs(Array1D(XReg))*Weight1D(XReg))
             OpEnWei = OpEnWei + Weight1D(XReg)
           else
             OpMissWei = OpMissWei + 1.0
           end if
         else
           OpMiss    = OpMiss    + 1.0
           
           if (Weight1D(XReg).NE.MissVal) then
             OpMissWei = OpMissWei + Weight1D(XReg)
           else
             OpMissWei = OpMissWei + 1.0
           end if
         end if
       end do
       
       if (OpEn.GT.0) then
         if ((OpMiss/(OpEn+OpMiss)).LT.(MissAccept/100.0)) then
           ArrayUM (XYear,XSeason) = OpTotUM / OpEn
           ArrayUA (XYear,XSeason) = OpTotUA / OpEn
         end if
       end if
       
       if (OpEnWei.GT.0) then
         if ((OpMissWei/(OpEnWei+OpMissWei)).LT.(MissAccept/100.0)) then
           ArrayWM (XYear,XSeason) = OpTotWM / OpEnWei
           ArrayWA (XYear,XSeason) = OpTotWA / OpEnWei
         end if
       end if
     end do
   end do
   
   SaveStem = trim(TextDir) // StatsText(XStat) // "-" // trim(TextTarg)
   if (XStat.NE.2) SaveStem = trim(SaveStem) // "." // trim(TextPatt)
   SaveStem = trim(SaveStem) // "." // trim(TextRegSet(XDumpRegSet))
   
   SaveTip = trim(CurrentExec%Name) // ".per"
   
   SaveFile = trim(SaveStem) // "-UM" // trim(SaveTip)
   call SavePer (SaveFile, YearAD, CallVariCode, AllData=ArrayUM, NoResponse=1)
   
   SaveFile = trim(SaveStem) // "-UA" // trim(SaveTip)
   call SavePer (SaveFile, YearAD, CallVariCode, AllData=ArrayUA, NoResponse=1)
   
   SaveFile = trim(SaveStem) // "-WM" // trim(SaveTip)
   call SavePer (SaveFile, YearAD, CallVariCode, AllData=ArrayWM, NoResponse=1)
   
   SaveFile = trim(SaveStem) // "-WA" // trim(SaveTip)
   call SavePer (SaveFile, YearAD, CallVariCode, AllData=ArrayWA, NoResponse=1)
   
  end if
 end do
 
 deallocate (MapIDLReg,RegSizes,RegNames,GloWeights,RegWeights,&
   		ArrayUM,ArrayUA,ArrayWM,ArrayWA,Array1D,Weight1D,Extract, stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpGlobe: Deallocation failure: Reg* #####"
end do

end subroutine DumpGlobe

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

subroutine PrepareRegions

call LoadRefReg (DumpRegSet(XDumpRegSet),MapIDLReg,RegSizes,RegNames,ExeN=ExeN,WyeN=WyeN)
RegN = size (RegNames)
 
call GetGloWeights (ExeN,WyeN,GloWeights)
 
allocate (RegWeights(RegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PrepareRegions: Allocation failure: RegWeights #####"
 
RegWeights = 0.0
do XExe = 1, ExeN
   do XWye = 1, WyeN
     if (MapIDLReg(XExe,XWye).NE.MissVal) then
       if (RegWeights(MapIDLReg(XExe,XWye)).NE.MissVal.AND.GloWeights(XExe,XWye).NE.MissVal) then
         RegWeights(MapIDLReg(XExe,XWye)) = RegWeights(MapIDLReg(XExe,XWye)) + GloWeights(XExe,XWye)
       else
         RegWeights(MapIDLReg(XExe,XWye)) = MissVal
       end if
     end if
   end do
end do
 
end subroutine PrepareRegions

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

subroutine MakeBoxVector

Extract = MissVal
if      (XStat.EQ.1) then
  do XBox = 1, BoxN
    Extract (XBox) = Est (XYear,XSeason,XBox)
  end do
else if (XStat.EQ.2) then
  do XBox = 1, BoxN
    Extract (XBox) = Mod (XYear,XSeason,XBox)
  end do
else if (XStat.EQ.3) then
  do XBox = 1, BoxN
    if (Est(XYear,XSeason,XBox).NE.MissVal.AND.Mod(XYear,XSeason,XBox).NE.MissVal) &
    			Extract (XBox) = Est(XYear,XSeason,XBox) - Mod(XYear,XSeason,XBox)
  end do
end if

end subroutine MakeBoxVector

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

subroutine MakeRegVector

Array1D = 0.0 ; Weight1D = 0.0
        
do XExe = 1, ExeN
  do XWye = 1, WyeN
             if (MapIDLReg(XExe,XWye).NE.MissVal) then
               if (Grid(XExe,XWye).NE.MissVal.AND.GloWeights(XExe,XWye).NE.MissVal) then
                if (Extract(Grid(XExe,XWye)).NE.MissVal) then
                  Array1D(MapIDLReg(XExe,XWye)) = Array1D(MapIDLReg(XExe,XWye)) + &
                  					(GloWeights(XExe,XWye)*Extract(Grid(XExe,XWye)))
                  Weight1D(MapIDLReg(XExe,XWye)) = Weight1D(MapIDLReg(XExe,XWye)) + GloWeights(XExe,XWye)
                end if
               end if
             end if            
  end do
end do
        
do XReg = 1, RegN
          if (Array1D(XReg).NE.MissVal) then
            if (Weight1D(XReg).GT.(RegWeights(XReg)*((100.0-MissAccept)/100.0))) then
          	Array1D(XReg) = Array1D(XReg) / Weight1D(XReg)
            else
          	Array1D(XReg) = MissVal
            end if
          end if
end do

end subroutine MakeRegVector

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

subroutine ResetCurrentExec

do						! resets CurrentExec so that it points at first Exec
      if      (associated(CurrentExec%Prev)) CurrentExec => CurrentExec%Prev
      if (.not.associated(CurrentExec%Prev)) exit
end do

end subroutine ResetCurrentExec

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

subroutine RemoveEstMod

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

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

end subroutine RemoveEstMod

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

subroutine RemoveExecs

call ResetCurrentExec

do
      if (associated(CurrentExec%Next)) then			! set stack to next Exec...
      	StackExec => CurrentExec%Next	
      else							! ...or if no next Exec, nullify
        nullify (StackExec)
      end if
      
      deallocate (CurrentExec, stat=AllocStat)			! deallocate current Exec
      if (AllocStat.NE.0) print*, "  > ##### ERROR: RemoveExecs: Deallocation failure: CurrentExec #####"
      
      if (     associated(StackExec)) CurrentExec => StackExec  ! set current to next Exec...
      if (.not.associated(StackExec)) exit			! ....or exit
end do

end subroutine RemoveExecs

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

subroutine Finish

deallocate (YearAD,DumpYear,DumpRegSet,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Finish: Deallocation failure #####"

call RemoveExecs

print*
close (99)

end subroutine Finish

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

end program ScalePattern