! seqgrim.f90
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
! 	-o ./../grim/seqgrim time.f90 filenames.f90 grimfiles.f90 annfiles.f90 
!   	initialmod.f90 sortmod.f90 plainperm.f90 regress.f90 linfiles.f90 
!   	glofiles.f90 saveperfiles.f90 aggfiles.f90 basicfun.f90 koeppen.f90 
!	./../grim/seqgrim.f90 2> /tyn1/tim/scratch/stderr.txt
! written by Tim Mitchell on 30.03.01
! last modified on 10.01.02
! it is designed to be run using Comms defined using a specification file, 
!	written in SetUpSeqGrim

program SeqGrim

use Time
use FileNames
use GrimFiles
use AnnFiles
use InitialMod
use SortMod
use PlainPerm
use Regress
use LinFiles
use GloFiles
use SavePerFiles
use AggFiles
use BasicFun
use Koeppen

implicit none
							! #######################
character (len=80), parameter :: SpecFile="./../../../code/f90/grim/_ref/130.ops"

real, pointer, dimension (:,:,:) :: OperA, OperB, OperC, OperD, OperK
real, pointer, dimension (:,:,:) :: Grim1, Grim2, Grim3, Grim4, Grim5, Grim6
real, pointer, dimension (:,:,:) :: GrimXtra
real, pointer, dimension (:,:)	 :: GridWeights, WeightsXtra, LinData
real, pointer, dimension (:,:)	 :: MonthData,SeasonData,MonthStDev,SeasonStDev
real, pointer, dimension (:,:)	 :: MonthTemp,MonthPrec
real, pointer, dimension (:)	 :: CommKay, VecA, VecB, VecC, VecD, GloSlice
real, pointer, dimension (:)	 :: RegBoreal,RegWeights,GloWeights,AnnualData,AnnualStDev
real, dimension (12)		 :: BoxTemp,BoxPrec
real, dimension (4)		 :: Bounds, FileBounds

integer, pointer, dimension (:,:) 		:: MapIDLReg	! shape of IDL mapping arrays
integer, pointer, dimension (:) 		:: RegSizes	! region sizes + reg-->raw
integer, pointer, dimension (:,:)		:: CommOp, Grid, FileGrid, MonLen
integer, pointer, dimension (:)			:: YearAD, FileYearAD

character (len=80), pointer, dimension (:,:)	:: CommFile, CommInfo
character (len=80), pointer, dimension (:)	:: FileLineNames,RegFilePaths
character (len=20), pointer, dimension (:) 	:: RegNames	!names of individual regions
character (len=09), pointer, dimension (:)	:: ColTitles

character (len= 2), dimension (12) :: MonthName = (/'01','02','03','04','05','06',&
						    '07','08','09','10','11','12'/)
character (len= 3), dimension (4) :: SeasonName = (/'MAM','JJA','SON','DJF'/)

real, parameter 		:: MissVal = -999.0
real, parameter 		:: MissAccept = 10.0
integer, parameter		:: SeasonN = 4

character (len=80), parameter 	:: Blank    = ""

logical :: Boreal,QIntercept

real :: OpTotSq, OpTot, OpEn, OpMiss, OpMean, OpThresh, OpValid, OpBelow, OpAbove
real :: FileFactor, LSRAye,LSRBee,LSRPea,LSRAre,LSREnn, NorthSouth,BoxLat

integer :: AllocStat, CheckGrid
integer :: YearN, MonthN, BoxN, ExeN, WyeN, ExecN, ArrayN, CommN, FileYearN, RegN
integer :: XYear, XMonth, XBox, XExe, XWye, XExec, XArray, XComm, XElement, XPerYear, XBound, XSeason
integer :: XFileYear, XMasterYear, XReg, XRegCheck, XLetter, XExeNear,XWyeNear
integer :: AStart,AEnd,BStart,BEnd,XAyear,XBYear
integer :: SuffixStart, FileSubBeg, InfoSubBeg, RegNameLen, KopNumb, QMeanSum
integer :: PerYear0, PerYear1, YearAD0, YearAD1, FileYear0, FileYear1, MasterYear0, MasterYear1
integer :: ThisMonth, ThisYear, NoZip, TotA,TotB,TotC,TotD, Operate, Year0,Year1

character (len=80) :: FileInfo,FilePath,FileVariable,GridFilePath,GloFile,GloInfo,GenericFile,SaveVariName
character (len=80) :: Stem,Tip,ExtraStr,GivenFile,InfoLine
character (len=20) :: RegName,LineFormat
character (len=40) :: RegTitle
character (len=30) :: KopName
character (len= 4) :: FileSuffix, SaveSuffix
character (len= 3) :: KopAcro

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

call Initialise
call Main
call Finalise

contains

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

subroutine Initialise

open  (99,file="./../../../scratch/log-seqgrim.dat",access="sequential", &
		form="formatted",status="replace",action="write")

nullify (OperA,OperB,OperC,OperD,OperK)
nullify (Grim1,Grim2,Grim3,Grim4,Grim5,Grim6)
nullify (GrimXtra,LinData,CommKay,VecA,VecB,VecC)
nullify (CommOp,Grid,FileGrid,YearAD,FileYearAD)
nullify (CommFile,CommInfo,FileLineNames)

print*
print*, "  > ***** SeqGrim: performs pre-determined sequence of grim ops *****"
print*
print "(2a)", "   > ", trim(SpecFile)
print* 

call LoadSpec

if (ArrayN.GE.1.AND.AllocStat.EQ.0) then
  allocate (Grim1 (YearN,MonthN,BoxN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Initialise: Allocation failure: Grim1 #####"
  
  if (ArrayN.GE.2.AND.AllocStat.EQ.0) then
    allocate (Grim2 (YearN,MonthN,BoxN), stat=AllocStat)
    if (AllocStat.NE.0) print*, "  > ##### ERROR: Initialise: Allocation failure: Grim2 #####"
  
    if (ArrayN.GE.3.AND.AllocStat.EQ.0) then
      allocate (Grim3 (YearN,MonthN,BoxN), stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: Initialise: Allocation failure: Grim3 #####"
  
      if (ArrayN.GE.4.AND.AllocStat.EQ.0) then  
        allocate (Grim4 (YearN,MonthN,BoxN), stat=AllocStat)
        if (AllocStat.NE.0) print*, "  > ##### ERROR: Initialise: Allocation failure: Grim4 #####"
  
        if (ArrayN.GE.5.AND.AllocStat.EQ.0) then  
          allocate (Grim5 (YearN,MonthN,BoxN), stat=AllocStat)
          if (AllocStat.NE.0) print*, "  > ##### ERROR: Initialise: Allocation failure: Grim5 #####"
  
          if (ArrayN.GE.6.AND.AllocStat.EQ.0) then  
            allocate (Grim6 (YearN,MonthN,BoxN), stat=AllocStat)
            if (AllocStat.NE.0) print*, "  > ##### ERROR: Initialise: Allocation failure: Grim6 #####"
          end if  
        end if  
      end if  
    end if  
  end if
end if

end subroutine Initialise

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

subroutine LoadSpec

write (99,"(a)"), "subroutine LoadSpec"
open  (1,file=SpecFile,status="old",access="sequential",form="formatted",action="read")

read (1,"(5i9)"), YearN, MonthN, BoxN, ExeN, WyeN
read (1,"(5i9)"), ExecN, ArrayN, CommN, YearAD0, YearAD1
read (1,"(4f9.2)"), (Bounds(XBound), XBound=1,4)

print*, "  > Basic specifications: "
print "(a13,i4,a11,i2,a11,i6,a10,i4,a9,i4)", "   > Years = ", YearN, ", Months = ", MonthN, &
			", Domain = ", BoxN, ", Longs = ", ExeN, ", Lats = ", WyeN
print "(a18,i2,a11,i1,a13,i2,a11,i4,a1,i4)", "   > Executions = ", ExecN, ", Arrays = ", ArrayN, &
			", Commands = ", CommN, ", Period = ", YearAD0, "-", YearAD1

allocate (CommOp   (CommN,5),     &
	  CommKay  (CommN),       &
	  CommFile (CommN,ExecN), &
	  CommInfo (CommN,ExecN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadSpec: Allocation failure: Comm* #####"

write (99,"(f10.2)") MissVal
write (99,"(a)") "SeqGrim : data found within comms file."
do XComm = 1, CommN
  read (1,"(5i9,f9.2)"),  (CommOp(XComm,XElement), XElement=1,5), CommKay(XComm)
  
  do XExec = 1, ExecN
    read (1,"(a80)"), CommFile(XComm,XExec)
    read (1,"(a80)"), CommInfo(XComm,XExec)
  end do
end do
write (99,"(a)") "SeqGrim : data found within comms file. END."

allocate (YearAD (YearN),     &
	  Grid   (ExeN,WyeN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadSpec: Allocation failure: G+YAD #####"

do XExe = 1, ExeN
  do XWye = 1, WyeN
    read (1,"(i10)"), Grid(XExe,XWye)
  end do
end do

close (1)

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

end subroutine LoadSpec

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

subroutine Main

do XExec = 1, ExecN
  print*
  call ResetAll
  call FindVariable
  
  do XComm = 1, CommN
    if      (CommOp(XComm,1).EQ.-1) then
        print "(a,i1)", "   > Resetting: ", CommOp(XComm,2)
        call ResetGrim
    else if (CommOp(XComm,1).EQ. 0) then
        print "(a,i1)", "   > Loading into: ", CommOp(XComm,2)
        call GrabGrim
    else if (CommOp(XComm,1).EQ. 1) then
        print "(a,i1)", "   > Loading single year into: ", CommOp(XComm,2)
        call GrabGrimYear
    else if (CommOp(XComm,1).EQ. 2) then
        print "(a,i1)", "   > Loading .lin into: ", CommOp(XComm,2)
	call GrabLin
    else if (CommOp(XComm,1).EQ. 3) then
	call GrabAnn
        print "(a,i2,a,i1,a,i3,2(a,i4))", "   > Loaded .ann (col ", CommOp(XComm,4), ") into: ", CommOp(XComm,2), &
        			"; VALID:", TotA, "; COMMON: ",YearAD(AStart),"-",YearAD(AEnd)
    else if (CommOp(XComm,1).EQ. 4) then
        print "(a,i1,2(a,i4))", "   > Loading into: ", CommOp(XComm,2), " period=",CommOp(XComm,4),&
        			"-", CommOp(XComm,5)
        call GrabGrim
    else if (CommOp(XComm,1).EQ.10) then
    	call OpClim
        print "(a,i1,a,i4,a1,i4,a,i1)", "   > Climatology: ", CommOp(XComm,2), " = ", CommOp(XComm,4), "-", &
        			CommOp(XComm,5), " mean from ", CommOp(XComm,3)
    else if (CommOp(XComm,1).EQ.11) then
    	call OpGloAv
        print "(a,i1,a,i1)", "   > Global mean: ", CommOp(XComm,2), " = average at each year,month from ", &
        			CommOp(XComm,3)
    else if (CommOp(XComm,1).EQ.12) then
    	call OpAnnAv
        print "(a,i1,a,i1)", "   > Annual mean: ", CommOp(XComm,2), " = average at each year,box from ", &
        			CommOp(XComm,3)
    else if (CommOp(XComm,1).EQ.13) then
    	call OpAnomalise
        if (CommOp(XComm,3).EQ.0) then
          print "(2(a,i1),2(a,i4),2(a,i9))", "   > Anomalise: ", CommOp(XComm,2), " = ", &
        			CommOp(XComm,2), " minus the mean of ", CommOp(XComm,4), &
        			"-", CommOp(XComm,5), " valid=", TotA, ", miss=", TotB
        else
          print "(2(a,i1),2(a,i4),2(a,i9))", "   > Anomalise: ", CommOp(XComm,2), " = ", &
        			CommOp(XComm,2), " minus the mean of ", CommOp(XComm,4), &
        			"-", CommOp(XComm,5), " valid=", TotA, ", igno=", TotC
        end if
    else if (CommOp(XComm,1).EQ.14) then
        if (CommOp(XComm,5).EQ.1) ExtraStr=" with missvals replaced"
        if (CommOp(XComm,5).EQ.2) ExtraStr=" with missvals retained"
        print "(a,i1,a,i1,a,i3,2a)", "   > Smooth: ", CommOp(XComm,2), " = ", CommOp(XComm,3), &
        			" smoothed at ", CommOp(XComm,4), " years", trim(ExtraStr)
    	call OpSmooth
    	print "(a,i10,a,i10)", "   >   VALID: ", TotA, " = ", TotB
    else if (CommOp(XComm,1).EQ.15) then
    	call OpConvert
        print "(a,i1,a,i1,a,i5,a,i5,a,i9)", "   > Convert: ", CommOp(XComm,2), " = ", CommOp(XComm,3), &
        	" but ", CommOp(XComm,4), "=>", CommOp(XComm,5), " ; changed: ", TotA
    else if (CommOp(XComm,1).EQ.16) then
    	call OpMask
        print "(a,i1,a,i5,a,i1,a,f8.3,a,i9)", "   > Where ", CommOp(XComm,2), "=", CommOp(XComm,4), &
        	", ", CommOp(XComm,3), "=", CommKay(XComm), " ; changed: ", TotA
    else if (CommOp(XComm,1).EQ.17) then
    	call OpMaskRel
        print "(a,i1,a,i2,3(a,i1),a,i9)", "   > Where ", CommOp(XComm,2), "[rel=", &
        	nint(CommKay(XComm)), "]", CommOp(XComm,5), &
        	", ", CommOp(XComm,3), "=", CommOp(XComm,4), " ; changed: ", TotA
    else if (CommOp(XComm,1).EQ.18) then
    	call OpMaskTransfer
        print "(a,i1,a,f8.3,a,i1,a,i5,a,i9)", "   > Where ", CommOp(XComm,2), "=", CommKay(XComm), &
        	", ", CommOp(XComm,3), "=", CommOp(XComm,4), " ; changed: ", TotA
    else if (CommOp(XComm,1).EQ.19) then
    	call OpMaskRelConstant
        print "(a,i1,a,i2,a,f8.3,a,i1,a,f8.3,a,i9)", "   > Where ", CommOp(XComm,2), "[rel=", &
        	CommOp(XComm,5), "]", CommKay(XComm), &
        	", ", CommOp(XComm,3), "=", CommKay(XComm), " ; changed: ", TotA
    else if (CommOp(XComm,1).EQ.20) then
        if (CommOp(XComm,3).EQ.MissVal) then        
          if (CommOp(XComm,4).NE.CommOp(XComm,5)) then
            print "(a,i1,a,i1,a,i1)", "   > Regressing (y=ax): x=", CommOp(XComm,5), &
        			", y=", CommOp(XComm,4), ": a=", CommOp(XComm,2)
    	  else
            print "(a,i1,a,i1)", "   > Regressing (y=at): y=", CommOp(XComm,4), ": a=", CommOp(XComm,2)
    	  end if
    	else
          if (CommOp(XComm,4).NE.CommOp(XComm,5)) then
            print "(a,i1,a,i1,a,i1,a,i1)", "   > Regressing (y=ax+b): x=", CommOp(XComm,5), &
        			", y=", CommOp(XComm,4), ": a=", CommOp(XComm,2), ": b=", CommOp(XComm,3)
    	  else
            print "(a,i1,a,i1,a,i1)", "   > Regressing (y=at+b): y=", CommOp(XComm,4), &
            			": a=", CommOp(XComm,2), ": b=", CommOp(XComm,3)
    	  end if
    	end if
    	
    	call OpLSR
    	
    	print "(a,3i10)", "   >   VALID (x,y,a): ", TotC,TotD,TotA
    else if (CommOp(XComm,1).EQ.21) then
        print "(a,i1)", "   > Putting month lengths into: ", CommOp(XComm,2)
	call GrabMonLen
    else if (CommOp(XComm,1).EQ.22) then
        print "(a,i1,a,i1,a,i5,a,i5)", "   > Placing ", CommOp(XComm,3), " into ", CommOp(XComm,2), &
        			" with range truncated to ", CommOp(XComm,4), " - ", CommOp(XComm,5) 
	call OpTruncate
    else if (CommOp(XComm,1).EQ.30) then
    	call OpGrimKay
        print "(3(a,i1),2(a,i10),a)", "   > Operating: ", CommOp(XComm,2), " = f(", CommOp(XComm,3), &
        			"): op ", CommOp(XComm,4), "; VALID: ",TotA,"=f(",TotB,")"
    else if (CommOp(XComm,1).EQ.31) then
    	call OpGrimGrim
        print "(4(a,i1),3(a,i10),a)", "   > Operating: ", CommOp(XComm,2), " = f(", CommOp(XComm,3), &
        			",", CommOp(XComm,5), "): op ", CommOp(XComm,4), &
        			"; VALID: ",TotA,"=f(",TotB,",",TotC,")"
    else if (CommOp(XComm,1).EQ.32) then
    	call OpSetToKay
        print "(a,i1,a,f7.2)", "   > Setting: ", CommOp(XComm,2), "=", CommKay(XComm)
    else if (CommOp(XComm,1).EQ.33) then
    	call OpGrimGrimExt
        print "(4(a,i1),3(a,i10),a)", "   > Operating: ", CommOp(XComm,2), " = f(", CommOp(XComm,3), &
        			",", 0, "): op ", CommOp(XComm,4), &
        			"; VALID: ",TotA,"=f(",TotB,",",TotC,")"
    else if (CommOp(XComm,1).EQ.43) then
        print "(a,i1)", "   > Saving from: ", CommOp(XComm,2)
	call DumpGrim
    else if (CommOp(XComm,1).EQ.44) then
        print "(a,i4,a1,i4,a,i1)", "   > Saving period (", CommOp(XComm,4), "-", CommOp(XComm,5), &
        			") from: ", CommOp(XComm,2)
	call DumpGrimPer
    else if (CommOp(XComm,1).EQ.45) then
        if (CommOp(XComm,4).EQ.1) print "(a,i4,a,i1,a)", "   > Saving year ", CommOp(XComm,3), &
        		" to .glo from ", CommOp(XComm,2), ": annual moment"
        if (CommOp(XComm,4).EQ.2) print "(a,i4,a,i1,a)", "   > Saving year ", CommOp(XComm,3), &
        		" to .glo from ", CommOp(XComm,2), ": seasonal moments"
        if (CommOp(XComm,4).EQ.3) print "(a,i4,a,i1,a)", "   > Saving year ", CommOp(XComm,3), &
        		" to .glo from ", CommOp(XComm,2), ": monthly values"
	call DumpGlo
    else if (CommOp(XComm,1).EQ.46) then
        RegName = ""					! all this is to get RegName=blank and SaveVariName
        GenericFile  = CommFile(XComm,XExec)
	SuffixStart  = index(GenericFile,".per")
	GenericFile  = GenericFile(1:(SuffixStart-1))
	SaveSuffix   = GetFileSuffix (GenericFile)
	call CheckVariSuffix (SaveSuffix,SaveVariName,FileFactor)

        print "(a,i6,a,i1)", "   > Saving box ", CommOp(XComm,3), " to .per from ", CommOp(XComm,2)
        call PointOperA		! don't alter without first understanding why this is here, from option 7
        call DumpDotPer        
	nullify (OperA)		! as above
    else if (CommOp(XComm,1).EQ.47) then
        print "(a,i1)", "   > Saving each region to .per from ", CommOp(XComm,2)
	call DumpRegPer
    else if (CommOp(XComm,1).EQ.48) then
	print "(a,2(i4,a),i1)", "   > Saving ", CommOp(XComm,3), "-", CommOp(XComm,5), &
			" to .agg from ", CommOp(XComm,2)
        call DumpAgg
    else if (CommOp(XComm,1).EQ.49) then
	print "(a,i1,a,f8.2,a)", "   > Counting ", CommOp(XComm,2), "=", CommKay(XComm), " to .glo"
        call CountOccurrences
    else if (CommOp(XComm,1).EQ.50) then
	print "(a,i1)", "   > Over-riding auto: setting QMeanSum = ", CommOp(XComm,2)
        QMeanSum = CommOp(XComm,2)
    else if (CommOp(XComm,1).EQ.51) then
        call InterpolateMissing
	print "(2(a,i1),2(a,i9))", "   > Interpol missing ", CommOp(XComm,2), "=", CommOp(XComm,3), &
			": done=", TotA, ", fail=", TotB
    else if (CommOp(XComm,1).EQ.52) then
        call MakeKoeppen
	print "(2(a,i1),2(a,i9))", "   > Koeppen from ", CommOp(XComm,2), "=temp ", CommOp(XComm,3), &
			"=prec ; done=", TotA, ", fail=", TotB
    end if
  end do
end do

end subroutine Main

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

subroutine ResetAll

if (associated(Grim1)) Grim1=MissVal
if (associated(Grim2)) Grim2=MissVal
if (associated(Grim3)) Grim3=MissVal
if (associated(Grim4)) Grim4=MissVal
if (associated(Grim5)) Grim5=MissVal
if (associated(Grim6)) Grim6=MissVal

end subroutine ResetAll

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

subroutine FindVariable

XComm = 0 ; FileSuffix = "" ; FileVariable = "" ; FileFactor = MissVal
do
  XComm = XComm + 1
  GivenFile = CommFile(XComm,XExec)
  
  if (GivenFile.NE."") then
  	SuffixStart = index(GivenFile,'.',.TRUE.)
  	if (GivenFile(SuffixStart:(SuffixStart+1)).EQ.".Z") then
  	  GivenFile(SuffixStart:(SuffixStart+1)) = "  "
  	  SuffixStart = index(GivenFile,'.',.TRUE.)
  	end if
  	if (GivenFile(SuffixStart:(SuffixStart+3)).EQ.".lin".OR. &
  		GivenFile(SuffixStart:(SuffixStart+3)).EQ.".ann") then
  	  GivenFile(SuffixStart:(SuffixStart+3)) = "    "
  	  SuffixStart = index(GivenFile,'.',.TRUE.)
  	end if
  	
  	if (SuffixStart.GT.0) then
  	  FilePath = CommFile(XComm,XExec)
  	  FileSuffix = FilePath(SuffixStart:(SuffixStart+3))  	
  	  call CheckVariSuffix (FileSuffix,FileVariable,FileFactor)
  	end if
  end if
  
  if (XComm.EQ.CommN.OR.FileVariable.NE."") exit
end do

if (FileVariable.EQ."") FileVariable = "unknown"
print "(2a)", "   > Variable: ", trim(FileVariable)

if (FileSuffix.EQ.".pre".OR.FileSuffix.EQ.".frs".OR.FileSuffix.EQ.".wet") then
  QMeanSum = 2
else
  QMeanSum = 0
end if

end subroutine FindVariable

!*******************************************************************************
! command = -1

subroutine ResetGrim

if      (CommOp(XComm,2).EQ.1) then
  Grim1 = MissVal
else if (CommOp(XComm,2).EQ.2) then
  Grim2 = MissVal
else if (CommOp(XComm,2).EQ.3) then
  Grim3 = MissVal
else if (CommOp(XComm,2).EQ.4) then
  Grim4 = MissVal
else if (CommOp(XComm,2).EQ.5) then
  Grim5 = MissVal
else if (CommOp(XComm,2).EQ.6) then
  Grim6 = MissVal
end if

end subroutine ResetGrim

!*******************************************************************************
! command=0 or =4
! if this fails under =4 at LoadGrim, it may be because the spec period is inappropriate

subroutine GrabGrim

call PointOperA
if      (CommOp(XComm,1).EQ.0) then
  call LoadGrim (OperA,FileGrid,FileYearAD,FileBounds,FileInfo,CommFile(XComm,XExec),"    ",FileSuffix,&
		MasterYearAD=YearAD)
else if (CommOp(XComm,1).EQ.4) then
  call LoadGrim (OperA,FileGrid,FileYearAD,FileBounds,FileInfo,CommFile(XComm,XExec),"    ",FileSuffix,&
		MasterYearAD=YearAD,OverRide=CommOp(XComm,4))
end if

call GridMismatch
nullify (OperA)

deallocate (FileGrid,FileYearAD,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabGrim: Deallocation failure: File* #####"

end subroutine GrabGrim

!*******************************************************************************
! command = 1 

subroutine GrabGrimYear

call PointOperA
call LoadGrim (GrimXtra,FileGrid,FileYearAD,FileBounds,FileInfo,CommFile(XComm,XExec),"    ",FileSuffix)

call GridMismatch
if (size(FileYearAD).NE.1) print*, "  > ##### ERROR: file period > 1 year #####"

deallocate (FileGrid,FileYearAD,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabGrimYear: Deallocation failure: File* #####"

do XMonth = 1, MonthN
  do XBox = 1, BoxN
    do XYear = 1, YearN
	OperA(XYear,XMonth,XBox) = GrimXtra(1,XMonth,XBox)
    end do
  end do
end do

nullify (OperA)

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

end subroutine GrabGrimYear

!*******************************************************************************
! command = 2

subroutine GrabLin

call PointOperA
call LoadLin (CommFile(XComm,XExec),FileYearAD,FileLineNames,LinData)

if (size(FileYearAD).NE.YearN)  print*, "  > ##### ERROR: loaded YearAD wrong length"
if      (CommOp(XComm,3).EQ.1) then
  if (FileYearAD(1).NE.YearAD(1)) print*, "  > ##### ERROR: loaded YearAD wrong years"
else if (CommOp(XComm,3).EQ.2) then
  ! no check, as the period is allowed to be anything - it is just ignored
end if

if (size(LinData,1).EQ.1) then
 do XYear = 1, YearN
  do XMonth = 1, MonthN 
    do XBox = 1, BoxN
      OperA(XYear,XMonth,XBox) = LinData(1,XYear)
    end do
  end do
 end do
else
  print*, "  > ##### ERROR: GrabLin: LinData size =/ 1 #####"
end if

nullify (OperA)

deallocate (FileYearAD,FileLineNames,LinData,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabLin: Deallocation failure #####"

end subroutine GrabLin

!*******************************************************************************
! command = 3

subroutine GrabAnn

call PointOperA
call LoadANN (CommFile(XComm,XExec),FileYearAD,ColTitles,LinData)

call CommonVecPer (YearAD,FileYearAD,AStart,AEnd,BStart,BEnd)

TotA = 0
if      (AStart.EQ.-999.OR.BStart.EQ.-999) then
 print*, "  > ##### ERROR: loaded YearAD wrong years #####"
else if (size(LinData,2).LT.CommOp(XComm,4)) then
 print*, "  > ##### ERROR: loaded array does not contain desired column #####"
else
 do XAYear = AStart,AEnd
  XBYear = BStart + XAYear - 1
  if (LinData(XBYear,CommOp(XComm,4)).NE.MissVal) TotA=TotA+1
  do XMonth = 1, MonthN 
    do XBox = 1, BoxN
      OperA(XAYear,XMonth,XBox) = LinData(XBYear,CommOp(XComm,4))
    end do
  end do
 end do
end if

nullify (OperA)
deallocate (FileYearAD,ColTitles,LinData,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabAnn: Deallocation failure #####"

end subroutine GrabAnn

!*******************************************************************************
! command = 10

subroutine OpClim

call PointOperA
call PointOperB

PerYear0 = 0 ; PerYear1 = 0 ; XYear = 0
do
   XYear = XYear + 1
   if (YearAD(XYear).EQ.CommOp(XComm,4)) PerYear0 = XYear
   if (YearAD(XYear).EQ.CommOp(XComm,5)) PerYear1 = XYear
   if (PerYear1.GT.0.OR.XYear.EQ.YearN) exit
end do

OpThresh = real(PerYear1-PerYear0+1) * (MissAccept/100.0)

if (PerYear1.GE.PerYear0) then
  do XBox = 1, BoxN
   do XMonth = 1, MonthN
    OpTot = 0.0 ; OpMiss = 0.0 ; OpEn = 0.0
   
    do XYear = PerYear0, PerYear1
     if (OperB(XYear,XMonth,XBox).NE.MissVal) then
     	OpTot  = OpTot  + OperB(XYear,XMonth,XBox)
     	OpEn   = OpEn   + 1.0
     else
     	OpMiss = OpMiss + 1.0
     end if
    end do
    
    if (OpEn.GT.0) then
      OpMean = OpTot / OpEn
    else
!      OpMean = 0.0
      OpMean = MissVal
    end if
    
!    if (OpMiss.LE.OpThresh.AND.OpEn.GT.0) then
!        OpMean = OpTot / OpEn
!    else
!        OpMean = MissVal
!    end if
   
    do XYear = 1, YearN
        OperA(XYear,XMonth,XBox) = OpMean
    end do
   end do
  end do
end if

nullify (OperA,OperB)

end subroutine OpClim

!*******************************************************************************
! command = 11

subroutine OpGloAv

call PointOperA
call PointOperB

do XYear = 1, YearN
 do XMonth = 1, MonthN
  OpTot = 0.0 ; OpMiss = 0.0 ; OpEn = 0.0

  do XBox = 1, BoxN   
    if (OperB(XYear,XMonth,XBox).NE.MissVal) then
     	OpTot  = OpTot  + OperB(XYear,XMonth,XBox)
     	OpEn   = OpEn   + 1.0
    else
     	OpMiss = OpMiss + 1.0
    end if   
  end do
  
  if (OpMiss.EQ.0.AND.OpEn.GT.0) then
        OpMean = OpTot / OpEn
  else
        OpMean = MissVal
  end if
   
  do XBox = 1, BoxN
        OperA(XYear,XMonth,XBox) = OpMean
  end do
 end do
end do

nullify (OperA,OperB)

end subroutine OpGloAv
 
!*******************************************************************************
! command = 12

subroutine OpAnnAv

call PointOperA
call PointOperB

do XYear = 1, YearN
 do XBox = 1, BoxN   
  OpTot = 0.0 ; OpMiss = 0.0 ; OpEn = 0.0

  do XMonth = 1, MonthN
    if (OperB(XYear,XMonth,XBox).NE.MissVal) then
     	OpTot  = OpTot  + OperB(XYear,XMonth,XBox)
     	OpEn   = OpEn   + 1.0
    else
     	OpMiss = OpMiss + 1.0
    end if   
  end do
  
  if (OpMiss.EQ.0.AND.OpEn.GT.0) then
        OpMean = OpTot / OpEn
  else
        OpMean = MissVal
  end if
   
  do XMonth = 1, MonthN
        OperA(XYear,XMonth,XBox) = OpMean
  end do
 end do
end do

nullify (OperA,OperB)

end subroutine OpAnnAv
 
!*******************************************************************************
! command = 13

subroutine OpAnomalise

call PointOperA

PerYear0 = 0 ; PerYear1 = 0 ; XYear = 0
do
   XYear = XYear + 1
   if (YearAD(XYear).EQ.CommOp(XComm,4)) PerYear0 = XYear
   if (YearAD(XYear).EQ.CommOp(XComm,5)) PerYear1 = XYear
   if (PerYear1.GT.0.OR.XYear.EQ.YearN) exit
end do

OpThresh = real(PerYear1-PerYear0+1) * (MissAccept/100.0)

TotA=0 ; TotB=0 ; TotC=0 
if (PerYear1.GT.PerYear0) then
  allocate (MonthData(MonthN,BoxN),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure: OpAnomalise #####"
  
  do XBox = 1, BoxN					! calc climatology
   do XMonth = 1, MonthN
    OpTot = 0.0 ; OpMiss = 0.0 ; OpEn = 0.0
   
    do XYear = PerYear0, PerYear1
     if (OperA(XYear,XMonth,XBox).NE.MissVal) then
     	OpTot  = OpTot  + OperA(XYear,XMonth,XBox)
     	OpEn   = OpEn   + 1.0
     else
     	OpMiss = OpMiss + 1.0
     end if
    end do
    
    if (OpMiss.LE.OpThresh) then
        OpMean = OpTot / OpEn
    else
        OpMean = MissVal
    end if
    
    MonthData(XMonth,XBox) = OpMean
   end do
  end do
  
  do XMonth = 1, MonthN
   do XBox = 1, BoxN
    if (MonthData(XMonth,XBox).NE.MissVal) then
      do XYear = 1, YearN
        if (OperA(XYear,XMonth,XBox).NE.MissVal) then
      	  OperA(XYear,XMonth,XBox) = OperA(XYear,XMonth,XBox) - MonthData(XMonth,XBox)
      	  TotA = TotA + 1
        end if
      end do
    else
     if (CommOp(XComm,3).EQ.1) then
      do XYear = 1, YearN
        TotC = TotC + 1
      end do
     else
      do XYear = 1, YearN
        OperA(XYear,XMonth,XBox) = MissVal
        TotB = TotB + 1
      end do
     end if
    end if
   end do
  end do
  
  deallocate (MonthData,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure: OpAnomalise #####"
end if

nullify (OperA)

end subroutine OpAnomalise

!*******************************************************************************
! command = 14

subroutine OpSmooth

call PointOperA
call PointOperB

allocate (VecA(YearN),VecB(YearN),VecC(YearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure: OpSmooth: Vec* #####"

TotA=0 ; TotB=0
do XMonth = 1, MonthN
  do XBox = 1, BoxN
    do XYear = 1, YearN
      VecB(XYear) = OperB(XYear,XMonth,XBox)
      if (VecB(XYear).NE.MissVal) TotB=TotB+1
    end do
    VecA=MissVal ; VecC=MissVal
    call GaussSmooth (YearN,CommOp(XComm,4),CommOp(XComm,5),VecB,VecA,VecC)
    do XYear = 1, YearN
      if (VecA(XYear).NE.MissVal) TotA=TotA+1
      OperA(XYear,XMonth,XBox) = VecA(XYear)
    end do
  end do
end do

deallocate (VecA,VecB,VecC, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure: OpSmooth: Vec* #####"
nullify (OperA,OperB)

end subroutine OpSmooth

!*******************************************************************************
! command = 15

subroutine OpConvert

call PointOperA
call PointOperB

TotA=0
do XMonth = 1, MonthN
  do XBox = 1, BoxN
    do XYear = 1, YearN
        if (OperB(XYear,XMonth,XBox).EQ.real(CommOp(XComm,4))) then
          OperA(XYear,XMonth,XBox) = real(CommOp(XComm,5))
          TotA=TotA+1
        else
          OperA(XYear,XMonth,XBox) = OperB(XYear,XMonth,XBox)
        end if
    end do
  end do
end do

nullify (OperA,OperB)

end subroutine OpConvert

!*******************************************************************************
! command = 16

subroutine OpMask

call PointOperA
call PointOperB

TotA=0
do XMonth = 1, MonthN
  do XBox = 1, BoxN
    do XYear = 1, YearN
      if (OperA(XYear,XMonth,XBox).EQ.real(CommOp(XComm,4))) then
        OperB(XYear,XMonth,XBox) = CommKay(XComm)
        TotA=TotA+1
      end if
    end do
  end do
end do

nullify (OperA,OperB)

end subroutine OpMask

!*******************************************************************************
! command = 17

subroutine OpMaskRel

call PointOperA
call PointOperB
call PointOperC
call PointOperD

TotA=0 ; Operate=nint(CommKay(XComm))
do XMonth = 1, MonthN
  do XBox = 1, BoxN
    do XYear = 1, YearN
      if (OperA(XYear,XMonth,XBox).NE.MissVal.AND.OperC(XYear,XMonth,XBox).NE.MissVal) then      
        if      (Operate.EQ.-3.AND.OperA(XYear,XMonth,XBox).NE.OperC(XYear,XMonth,XBox)) then
          OperB(XYear,XMonth,XBox)=OperD(XYear,XMonth,XBox) ; TotA=TotA+1
        else if (Operate.EQ.-2.AND.OperA(XYear,XMonth,XBox).LT.OperC(XYear,XMonth,XBox)) then
          OperB(XYear,XMonth,XBox)=OperD(XYear,XMonth,XBox) ; TotA=TotA+1
        else if (Operate.EQ.-1.AND.OperA(XYear,XMonth,XBox).LE.OperC(XYear,XMonth,XBox)) then
          OperB(XYear,XMonth,XBox)=OperD(XYear,XMonth,XBox) ; TotA=TotA+1
        else if (Operate.EQ. 0.AND.OperA(XYear,XMonth,XBox).EQ.OperC(XYear,XMonth,XBox)) then
          OperB(XYear,XMonth,XBox)=OperD(XYear,XMonth,XBox) ; TotA=TotA+1
        else if (Operate.EQ. 1.AND.OperA(XYear,XMonth,XBox).GE.OperC(XYear,XMonth,XBox)) then
          OperB(XYear,XMonth,XBox)=OperD(XYear,XMonth,XBox) ; TotA=TotA+1
        else if (Operate.EQ. 2.AND.OperA(XYear,XMonth,XBox).GT.OperC(XYear,XMonth,XBox)) then
          OperB(XYear,XMonth,XBox)=OperD(XYear,XMonth,XBox) ; TotA=TotA+1
        end if
      end if
    end do
  end do
end do

nullify (OperA,OperB,OperC,OperD)

end subroutine OpMaskRel

!*******************************************************************************
! command = 18

subroutine OpMaskTransfer

call PointOperA
call PointOperB
call PointOperD

TotA=0
do XBox=1, BoxN
    do XMonth = 1, MonthN
     do XYear = 1, YearN
       if (OperA(XYear,XMonth,XBox).EQ.CommKay(XComm)) then
          OperB(XYear,XMonth,XBox) = OperD(XYear,XMonth,XBox)
          TotA=TotA+1
       end if
     end do
    end do
end do

nullify (OperA,OperB,OperD)

end subroutine OpMaskTransfer

!*******************************************************************************
! command = 19

subroutine OpMaskRelConstant

call PointOperA
call PointOperB

TotA=0 ; Operate=CommOp(XComm,5)
do XBox = 1, BoxN
  do XMonth = 1, MonthN
    do XYear = 1, YearN
      if (OperA(XYear,XMonth,XBox).NE.MissVal.OR.CommOp(XComm,4).EQ.1) then      
        if      (Operate.EQ.-3.AND.OperA(XYear,XMonth,XBox).NE.CommKay(XComm)) then
          OperB(XYear,XMonth,XBox)=CommKay(XComm) ; TotA=TotA+1
        else if (Operate.EQ.-2.AND.OperA(XYear,XMonth,XBox).LT.CommKay(XComm)) then
          OperB(XYear,XMonth,XBox)=CommKay(XComm) ; TotA=TotA+1
        else if (Operate.EQ.-1.AND.OperA(XYear,XMonth,XBox).LE.CommKay(XComm)) then
          OperB(XYear,XMonth,XBox)=CommKay(XComm) ; TotA=TotA+1
        else if (Operate.EQ. 0.AND.OperA(XYear,XMonth,XBox).EQ.CommKay(XComm)) then
          OperB(XYear,XMonth,XBox)=CommKay(XComm) ; TotA=TotA+1
        else if (Operate.EQ. 1.AND.OperA(XYear,XMonth,XBox).GE.CommKay(XComm)) then
          OperB(XYear,XMonth,XBox)=CommKay(XComm) ; TotA=TotA+1
        else if (Operate.EQ. 2.AND.OperA(XYear,XMonth,XBox).GT.CommKay(XComm)) then
          OperB(XYear,XMonth,XBox)=CommKay(XComm) ; TotA=TotA+1
        end if
      end if
    end do
  end do
end do

nullify (OperA,OperB)

end subroutine OpMaskRelConstant

!*******************************************************************************
! command = 20

subroutine OpLSR

if      (CommOp(XComm,3).EQ.MissVal) then
  QIntercept=.FALSE.
else
  QIntercept=.TRUE.
end if

call PointOperA
call PointOperC
call PointOperD
if (QIntercept) call PointOperB
if (CommKay(XComm).NE.MissVal) call PointOperK

allocate (VecC(YearN),VecD(YearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Allocation failure: OpLSR: Vec* #####"

TotA=0 ; TotC=0 ; TotD=0
do XMonth = 1, MonthN
  do XBox = 1, BoxN
    do XYear = 1, YearN
      if (CommOp(XComm,4).EQ.CommOp(XComm,5)) then
      	VecC(XYear) = real (YearAD(XYear))
      else
      	VecC(XYear) = OperC(XYear,XMonth,XBox)
      end if
      
      VecD(XYear) = OperD(XYear,XMonth,XBox)
      
      if (VecC(XYear).NE.MissVal) TotC=TotC+1
      if (VecD(XYear).NE.MissVal) TotD=TotD+1
    end do
    
    if (QIntercept) then		! y=ax+b
!       call LinearLSRVec (VecC,VecD,LSRBee,LSRAye,LSRPea)
       call RegressVectors (VecC,VecD,LSRBee,LSRAye,LSRAre,LSREnn,5)
       
       if (LSRAye.NE.MissVal) TotA=TotA+YearN
       do XYear = 1, YearN
         OperA(XYear,XMonth,XBox) = LSRAye
	 OperB(XYear,XMonth,XBox) = LSRBee
!	 if (CommKay(XComm).NE.MissVal) OperK(XYear,XMonth,XBox) = LSRPea
	 if (CommKay(XComm).NE.MissVal) OperK(XYear,XMonth,XBox) = LSRAre
       end do
    else				! y=ax
       call LinearLSRZeroVec (VecC,VecD,LSRAye)
       if (LSRAye.NE.MissVal) TotA=TotA+YearN
       do XYear = 1, YearN
         OperA(XYear,XMonth,XBox) = LSRAye
       end do
    end if
  end do
end do

deallocate (VecC,VecD, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Deallocation failure: OpLSR: Vec* #####"
nullify (OperA,OperC,OperD)
if (CommOp(XComm,3).NE.MissVal) nullify (OperB)
  
end subroutine OpLSR

!*******************************************************************************
! command = 21

subroutine GrabMonLen

allocate (MonLen(YearN,MonthN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Alloc: GrabMonLen #####"

call GetMonthLengths (YearAD,MonLen)

call PointOperA

do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XBox = 1, BoxN
      OperA(XYear,XMonth,XBox) = real(MonLen(XYear,XMonth))
    end do
  end do
end do

nullify (OperA)

deallocate (MonLen, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabMonLen: dealloc #####"

end subroutine GrabMonLen

!*******************************************************************************
! command = 22

subroutine OpTruncate

call PointOperA
call PointOperB

OpBelow = 0.0 ; OpAbove = 0.0 ; OpMiss = 0.0

do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XBox = 1, BoxN
      if (OperB(XYear,XMonth,XBox).NE.MissVal) then
       if (CommOp(XComm,4).NE.MissVal) then
        if (OperB(XYear,XMonth,XBox).GE.CommOp(XComm,4)) then
            OperA(XYear,XMonth,XBox) = OperB(XYear,XMonth,XBox)
        else
            OperA(XYear,XMonth,XBox) = CommOp(XComm,4)
            OpBelow = OpBelow + 1.0
        end if
       end if
      
       if (CommOp(XComm,5).NE.MissVal) then
        if (OperB(XYear,XMonth,XBox).LE.CommOp(XComm,5)) then
            OperA(XYear,XMonth,XBox) = OperB(XYear,XMonth,XBox)
        else
            OperA(XYear,XMonth,XBox) = CommOp(XComm,5)
            OpAbove = OpAbove + 1.0
        end if
       end if
      else
       OpMiss = OpMiss + 1.0
      end if
    end do
  end do
end do

nullify (OperA,OperB)

OpTot   = YearN * MonthN * BoxN
OpValid = OpTot - OpMiss

if (OpTot.GT.0.AND.OpValid.GT.0) then
  OpMiss  = 100.0 * OpMiss  / OpTot
  OpBelow = 100.0 * OpBelow / OpValid
  OpAbove = 100.0 * OpAbove / OpValid
else
  OpMiss = MissVal ; OpBelow = MissVal ; OpAbove = MissVal
end if

print "(a,f14.1,f8.4)", "   > Total data, and % missing: ", OpTot, OpMiss
print "(a,2f8.4)", "   > Of valid data, % truncated from < and > range: ", OpBelow, OpAbove

end subroutine OpTruncate

!*******************************************************************************
! command = 30

subroutine OpGrimKay

call PointOperA
call PointOperB

TotA=0 ; TotB=0 ; TotC=0
do XYear = 1, YearN
  do XMonth = 1, MonthN
   do XBox = 1, BoxN
      if (OperB(XYear,XMonth,XBox).NE.MissVal) TotB=TotB+1
      OperA(XYear,XMonth,XBox) = OpAwithB (OperB(XYear,XMonth,XBox),CommKay(XComm),CommOp(XComm,4))
      if (OperA(XYear,XMonth,XBox).NE.MissVal) TotA=TotA+1
   end do
  end do
end do

nullify (OperA,OperB)

end subroutine OpGrimKay

!*******************************************************************************
! command = 31

subroutine OpGrimGrim

call PointOperA
call PointOperB
call PointOperC

TotA=0 ; TotB=0 ; TotC=0
do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XBox = 1, BoxN
      if (OperB(XYear,XMonth,XBox).NE.MissVal) TotB=TotB+1
      if (OperC(XYear,XMonth,XBox).NE.MissVal) TotC=TotC+1
      OperA(XYear,XMonth,XBox) = OpAwithB (OperB(XYear,XMonth,XBox),OperC(XYear,XMonth,XBox),CommOp(XComm,4))
      if (OperA(XYear,XMonth,XBox).NE.MissVal) TotA=TotA+1
    end do
  end do
end do

nullify (OperA,OperB,OperC)

end subroutine OpGrimGrim

!*******************************************************************************
! command = 32

subroutine OpSetToKay

call PointOperA

do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XBox = 1, BoxN
      OperA(XYear,XMonth,XBox) = CommKay(XComm)
    end do
  end do
end do

nullify (OperA)

end subroutine OpSetToKay

!*******************************************************************************
! command = 33

subroutine OpGrimGrimExt

call PointOperA
call PointOperB

call LoadGrim (GrimXtra,FileGrid,FileYearAD,FileBounds,FileInfo,CommFile(XComm,XExec),"    ",FileSuffix)

call GridMismatch
if (size(FileYearAD).NE.1) print*, "  > ##### ERROR: file period > 1 year #####"

deallocate (FileGrid,FileYearAD,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: OpGrimGrimExt: Deallocation failure: File* #####"

TotA=0 ; TotB=0 ; TotC=0
do XYear = 1, YearN
  do XMonth = 1, MonthN
    do XBox = 1, BoxN
      if (OperB(XYear,XMonth,XBox).NE.MissVal) TotB=TotB+1
      if (GrimXtra(1,XMonth,XBox).NE.MissVal)  TotC=TotC+1
      OperA(XYear,XMonth,XBox) = OpAwithB (OperB(XYear,XMonth,XBox),GrimXtra(1,XMonth,XBox),CommOp(XComm,4))
      if (OperA(XYear,XMonth,XBox).NE.MissVal) TotA=TotA+1
    end do
  end do
end do

nullify (OperA,OperB)

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

end subroutine OpGrimGrimExt

!*******************************************************************************
! command = 43

subroutine DumpGrim

call PointOperA
call SaveGrim (OperA,Grid,YearAD,Bounds,CommInfo(XComm,XExec),CommFile(XComm,XExec),"    ",FileSuffix,NoZip=1)
nullify (OperA)

end subroutine DumpGrim

!*******************************************************************************
! command = 44

subroutine DumpGrimPer

call PointOperA

call SaveGrim (OperA,Grid,YearAD,Bounds,CommInfo(XComm,XExec),CommFile(XComm,XExec),"    ",&
		FileSuffix,NoZip=1,SaveYearAD0=CommOp(XComm,4),SaveYearAD1=CommOp(XComm,5))

nullify (OperA)

end subroutine DumpGrimPer

!*******************************************************************************
! command = 45

subroutine DumpGlo

GridFilePath = GetGloRef (ExeN,WyeN)
call PointOperA
  
XYear = CommOp(XComm,3) - YearAD(1) + 1
  
allocate (GloSlice(BoxN),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpGlo: Allocation failure."

if      (CommOp(XComm,4).EQ.1) then				! save annual mean
    GloSlice = MissVal
    
    do XBox = 1, BoxN
      OpTot = 0.0 ; OpEn = 0.0
      
      do XMonth = 1, MonthN
        if (OperA(XYear,XMonth,XBox).NE.MissVal) then
          OpTot = OpTot + OperA(XYear,XMonth,XBox)
          OpEn  = OpEn  + 1.0
        end if
      end do
      
      if (OpEn.EQ.MonthN) then
      	if (QMeanSum.EQ.0) GloSlice (XBox) = OpTot / OpEn
      	if (QMeanSum.EQ.2) GloSlice (XBox) = OpTot 
      end if
    end do
    
    call SaveGlo (ExeN,WyeN,BoxN,GridFilePath,CommFile(XComm,XExec),CommInfo(XComm,XExec),GloSlice,Grid)
        
else if (CommOp(XComm,4).EQ.2) then				! save seasonal means
    do XSeason = 1, 4						! note that XSeason=1 for MAM
      FileSubBeg = index(CommFile(XComm,XExec),SeasonName(1))
      if (FileSubBeg.GT.0) then
        GloFile = CommFile(XComm,XExec)
        GloFile (FileSubBeg:(FileSubBeg+2)) = SeasonName(XSeason)
      else
        GloFile = ""
      end if
      
      InfoSubBeg = index(CommInfo(XComm,XExec),SeasonName(1))
      if (InfoSubBeg.GT.0) then
        GloInfo = CommInfo(XComm,XExec)
        GloInfo (InfoSubBeg:(InfoSubBeg+2)) = SeasonName(XSeason)
      else
        GloInfo = ""
      end if
      
      if (GloFile.NE."".AND.GloInfo.NE."") then
        GloSlice = MissVal
        
        do XBox = 1, BoxN
          OpTot = 0.0 ; OpEn = 0.0
      
          do ThisMonth = (XSeason*3), (2+(XSeason*3))
            XMonth = ThisMonth
            if (XMonth.GT.12) XMonth = XMonth - 12
            
            if (OperA(XYear,XMonth,XBox).NE.MissVal) then
              OpTot = OpTot + OperA(XYear,XMonth,XBox)
              OpEn  = OpEn  + 1.0
            end if
          end do          
          
          if (OpEn.EQ.3) then
            if (QMeanSum.EQ.0) GloSlice (XBox) = OpTot / OpEn
      	    if (QMeanSum.EQ.2) GloSlice (XBox) = OpTot 
          end if
        end do
        
        call SaveGlo (ExeN,WyeN,BoxN,GridFilePath,GloFile,GloInfo,GloSlice,Grid)
      end if          
    end do
else if (CommOp(XComm,4).EQ.3) then				! save monthly values
    do XMonth = 1, MonthN
      FileSubBeg = index(CommFile(XComm,XExec),MonthName(1),.true.)
      if (FileSubBeg.GT.0) then
        GloFile = CommFile(XComm,XExec)
        GloFile (FileSubBeg:(FileSubBeg+1)) = MonthName(XMonth)
      else
        GloFile = ""
      end if
      
      InfoSubBeg = index(CommInfo(XComm,XExec),MonthName(1),.true.)
      if (InfoSubBeg.GT.0) then
        GloInfo = CommInfo(XComm,XExec)
        GloInfo (InfoSubBeg:(InfoSubBeg+1)) = MonthName(XMonth)
      else
        GloInfo = ""
      end if
      
      if (GloFile.NE."".AND.GloInfo.NE."") then
        do XBox = 1, BoxN
          GloSlice (XBox) = OperA(XYear,XMonth,XBox)
        end do
        
        call SaveGlo (ExeN,WyeN,BoxN,GridFilePath,GloFile,GloInfo,GloSlice,Grid)
      end if
    end do
end if          
  
deallocate (GloSlice,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpGlo: Deallocation failure."

nullify (OperA)

end subroutine DumpGlo

!*******************************************************************************
! command = 46 ; also called many times from command = 47

subroutine DumpDotPer

allocate (MonthData (YearN,MonthN),  &
	  SeasonData(YearN,SeasonN), &
	  AnnualData(YearN),         stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpDotPer: Allocation failure #####"
MonthData = MissVal ; SeasonData = MissVal ; AnnualData = MissVal

do XYear = 1, YearN				! fill MonthData
  do XMonth = 1, MonthN
    MonthData(XYear,XMonth) = OperA(XYear,XMonth,CommOp(XComm,3))
  end do
end do

do ThisYear = 1, YearN				! fill SeasonData
  do XSeason = 1, SeasonN  			! note that XSeason=1 for MAM
   OpTot = 0.0 ; OpEn = 0.0
   do ThisMonth = (XSeason*3), (2+(XSeason*3))
    XMonth = ThisMonth
    XYear  = ThisYear
    if (XMonth.GT.12) then
    	XMonth = XMonth - 12
!        XYear = XYear + 1	! command removed to enable single year dump
    end if
    
    if (XYear.LE.YearN) then
      if (OperA(XYear,XMonth,CommOp(XComm,3)).NE.MissVal) then
        OpTot = OpTot + OperA(XYear,XMonth,CommOp(XComm,3))
        OpEn  = OpEn  + 1
      end if
    end if
   end do
   if (OpEn.EQ.3) then
    if (QMeanSum.EQ.0) SeasonData(ThisYear,XSeason) = OpTot/OpEn
    if (QMeanSum.EQ.2) SeasonData(ThisYear,XSeason) = OpTot
   end if
  end do
end do

do XYear = 1, YearN				! fill AnnualData
  OpTot = 0.0 ; OpEn = 0.0
  do XMonth = 1, MonthN
    if (OperA(XYear,XMonth,CommOp(XComm,3)).NE.MissVal) then
      OpTot = OpTot + OperA(XYear,XMonth,CommOp(XComm,3))
      OpEn  = OpEn  + 1
    end if
  end do
  if (OpEn.EQ.MonthN) then
    if (QMeanSum.EQ.0) AnnualData(XYear) = OpTot/OpEn
    if (QMeanSum.EQ.2) AnnualData(XYear) = OpTot
  end if
end do

LineFormat='(i5,12f8.2,5f8.2)'
call SavePER (CommFile(XComm,XExec),YearAD,QMeanSum,Monthly=MonthData,Seasonal=SeasonData,&
		Annual=AnnualData,CallLineFormat=LineFormat,NoResponse=1)
!		ScenHeads=1,CallPoliUnitName=RegName,CallVariableName=SaveVariName)

deallocate (MonthData,SeasonData,AnnualData,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpDotPer: Deallocation failure #####"

end subroutine DumpDotPer

!*******************************************************************************
! command = 47

subroutine DumpRegPer

call PointOperA			! before we load up into regions, we point OperA towards the correct grim

call RegSelect (nint(MissVal), ExeN, WyeN, (ExeN*WyeN), MapIDLReg, RegSizes, RegNames, &
			RegTitle, RegN, CommInfo(XComm,1))
call GetWeightArrays

allocate (GrimXtra(YearN,MonthN,RegN), &
	  RegFilePaths(RegN),          stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpRegPer: Allocation failure #####"
GrimXtra = 0.0 				! careful - do not inintialise to MissVal

do XExe = 1, ExeN					! iterate around grid
  do XWye = 1, WyeN					! if the grid cell belongs to a region
    if (MapIDLReg(XExe,XWye).NE.MissVal) then		! then check if the grid cell is within the grim domain
      if (Grid(XExe,XWye).NE.MissVal.AND.RegWeights(MapIDLReg(XExe,XWye)).NE.MissVal) then		
        do XYear = 1, YearN				! iterate by year, month
          do XMonth = 1, MonthN
            if (OperA(XYear,XMonth,Grid(XExe,XWye)).NE.MissVal) then		! check datum = non-missing
              if (GrimXtra(XYear,XMonth,MapIDLReg(XExe,XWye)).NE.MissVal) then	! check region = OK thus far
                GrimXtra(XYear,XMonth,MapIDLReg(XExe,XWye)) = GrimXtra(XYear,XMonth,MapIDLReg(XExe,XWye)) &
          							+ (OperA(XYear,XMonth,Grid(XExe,XWye)) &
          							* GridWeights(XExe,XWye))
              end if
            else						! datum=missing, so make region missing
              GrimXtra(XYear,XMonth,MapIDLReg(XExe,XWye)) = MissVal
            end if
          end do
        end do
      else						
        do XYear = 1, YearN
          do XMonth = 1, MonthN
            GrimXtra(XYear,XMonth,MapIDLReg(XExe,XWye)) = MissVal  ! cell=missing, so make region missing
          end do
        end do
      end if
    end if
  end do
end do

do XReg = 1, RegN					! convert from sums to means
  do XYear = 1, YearN
    do XMonth = 1, MonthN
      if (RegWeights(XReg).GT.0.0) then
        if (GrimXtra(XYear,XMonth,XReg).NE.MissVal) then
          GrimXtra(XYear,XMonth,XReg) = GrimXtra(XYear,XMonth,XReg) / RegWeights(XReg)
        else
          GrimXtra(XYear,XMonth,XReg) = MissVal
        end if
      else
        GrimXtra(XYear,XMonth,XReg) = MissVal
      end if
    end do
  end do    
end do

nullify (OperA) ; OperA => GrimXtra	! we now point OperA towards the region array, ready for call to opt 6

GenericFile  = CommFile(XComm,XExec)			! generate region filepaths with correct directory
SuffixStart  = index(GenericFile,".per")
GenericFile  = GenericFile(1:(SuffixStart-1))
SaveSuffix   = GetFileSuffix (GenericFile)
SuffixStart  = index(GenericFile,SaveSuffix)

call CheckVariSuffix (SaveSuffix,SaveVariName,FileFactor)
Stem         = GenericFile(1:(SuffixStart-1)) // "."
Tip          = trim(SaveSuffix) // ".per"
call MakeBatch (Stem,Tip,RegNames,RegFilePaths)

do XReg = 1, RegN
  CommOp(XComm,3) = XReg
  CommFile(XComm,XExec) = RegFilePaths(XReg)
  RegName = RegNames(XReg)
  call DumpDotPer
end do

nullify (OperA)

deallocate (MapIDLReg, RegSizes, RegNames, GridWeights, RegWeights, RegBoreal, &
		GrimXtra, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpRegPer: Deallocation failure #####"

end subroutine DumpRegPer

!*******************************************************************************
! command = 48 - calcs mean and stdev and dumps separately

subroutine DumpAgg

call PointOperA		! before we load up into regions, we point OperA towards the correct grim

call RegSelect (nint(MissVal), ExeN, WyeN, (ExeN*WyeN), MapIDLReg, RegSizes, RegNames, &
			RegTitle, RegN, CommInfo(XComm,1))
call GetWeightArrays

Year0=CommOp(XComm,3)-YearAD(1)+1 ; Year1=CommOp(XComm,5)-YearAD(1)+1  
FileYearN=Year1-Year0+1		! 3,5 is deliberate

allocate (GrimXtra(FileYearN,17,RegN),   stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpAgg: Allocation failure: MonthData #####"
GrimXtra=0

print*, "  > Aggregating into regions..." 
do XYear = Year0,Year1
  XFileYear=XYear-Year0+1
  do XMonth = 1, MonthN
   do XExe=1,ExeN
    do XWye=1,WyeN
      XReg=MapIDLReg(XExe,XWye)
!      write (99,"(7i6)"), XYear,XFileYear,XMonth,XExe,XWye,XReg,Grid(XExe,XWye)
      if (XReg.NE.MissVal) then
!       write (99,"(2f10.2)"), GrimXtra(XFileYear,XMonth,XReg),OperA(XYear,XMonth,Grid(XExe,XWye))
       if (GrimXtra(XFileYear,XMonth,XReg).NE.MissVal) then
        if (Grid(XExe,XWye).NE.MissVal) then
          if (OperA(XYear,XMonth,Grid(XExe,XWye)).NE.MissVal) then
              GrimXtra(XFileYear,XMonth,XReg) = GrimXtra(XFileYear,XMonth,XReg) + &
              	(OperA(XYear,XMonth,Grid(XExe,XWye)) * GridWeights(XExe,XWye))
          else
              GrimXtra(XFileYear,XMonth,XReg) = MissVal
          end if
        end if
       end if		
      end if
    end do
   end do
  end do
end do

print*, "  > Removing aggregated weights..." 
do XReg=1,RegN
  if (RegWeights(XReg).GT.0) then
   do XfileYear=1,FileYearN
    do XMonth=1,MonthN
      if (GrimXtra(XFileYear,XMonth,XReg).NE.MissVal) &
      	GrimXtra(XFileYear,XMonth,XReg) = GrimXtra(XFileYear,XMonth,XReg) / &
      		RegWeights(XReg)
    end do
   end do
  end if
end do

print*, "  > Seasonalising regions..." 
do XReg=1,RegN		! fill GrimXtra with regional values for seasons/ years
  do XFileYear=1,FileYearN
    OpTot=0 ; OpEn=0
    do XMonth = 1, MonthN
      if (GrimXtra(XFileYear,XMonth,XReg).NE.MissVal) then
        OpTot=OpTot+GrimXtra(XFileYear,XMonth,XReg) ; opEn=OpEn+1
      end if
    end do
    if (OpEn.Eq.12) then
      if (QMeanSum.EQ.0) GrimXtra(XFileYear,17,XReg)=OpTot/OpEn
      if (QMeanSum.EQ.2) GrimXtra(XFileYear,17,XReg)=OpTot
    end if
    
    do XSeason = 1, SeasonN  			! note that XSeason=1 for MAM
      OpTot = 0.0 ; OpEn = 0.0
   
      do ThisMonth = (XSeason*3), (2+(XSeason*3))
        XMonth = ThisMonth ; XYear = XFileYear
        if (XMonth.GT.12) then		
          XMonth = XMonth - 12 
!          XYear = XYear + 1	! command removed to enable single year dump
        end if   
        if (XYear.LE.FileYearN) then
          if (GrimXtra(XYear,XMonth,XReg).NE.MissVal) then
            OpTot = OpTot + GrimXtra(XYear,XMonth,XReg)
            OpEn  = OpEn  + 1
          end if
        end if
      end do
      
      if (OpEn.Eq.3) then
        if (QMeanSum.EQ.0) GrimXtra(XFileYear,(12+XSeason),XReg)=OpTot/OpEn
        if (QMeanSum.EQ.2) GrimXtra(XFileYear,(12+XSeason),XReg)=OpTot
      end if
    end do
  end do
end do

nullify (OperA)

deallocate (MapIDLReg, RegSizes, GridWeights, RegWeights, RegBoreal, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpAgg: Deallocation failure: v1 #####"

allocate (SeasonData(RegN,SeasonN),MonthData(RegN,MonthN), &
	  AnnualData(RegN),SeasonStDev(RegN,SeasonN),MonthStDev(RegN,MonthN), &
	  AnnualStDev(RegN),         stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: NewDumpAgg: Allocation failure: SeaAnnStDev #####"
SeasonData=MissVal ; MonthData=MissVal ; AnnualData=MissVal
SeasonStDev=MissVal ; MonthStDev=MissVal ; AnnualStDev=MissVal

print*, "  > Calculating mean and stdev..." 
do XReg = 1, RegN			! calc mean and stdev
  do XMonth = 1, 17
    OpTot=0 ; OpTotSq=0 ; OpEn=0
    do XFileYear=1,FileYearN
      if (GrimXtra(XFileYear,XMonth,XReg).NE.MissVal) then
        OpTot=OpTot+GrimXtra(XFileYear,XMonth,XReg)
        OpTotSq=OpTotSq+(GrimXtra(XFileYear,XMonth,XReg)**2)
        OpEn=OpEn+1
      end if
    end do
    if (OpEn.GT.0) then
      if (XMonth.LE.12) then
        MonthData(XReg,XMonth)=OpTot/OpEn
        MonthStDev(XReg,XMonth)=(OpTotSq/OpEn)-(MonthData(XReg,XMonth)**2)
        if (MonthStDev(XReg,XMonth).GT.0) then
          MonthStDev(XReg,XMonth)=sqrt(MonthStDev(XReg,XMonth))
        else
          MonthStDev(XReg,XMonth)=0
        end if
      else if (XMonth.LE.16) then
        SeasonData(XReg,(XMonth-12))=OpTot/OpEn
        SeasonStDev(XReg,(XMonth-12))=(OpTotSq/OpEn)-(SeasonData(XReg,(XMonth-12))**2)
        if (SeasonStDev(XReg,(XMonth-12)).GT.0) then
          SeasonStDev(XReg,(XMonth-12))=sqrt(SeasonStDev(XReg,(XMonth-12)))
        else
          SeasonStDev(XReg,(XMonth-12))=0
        end if
      else
        AnnualData(XReg)=OpTot/OpEn
        AnnualStDev(XReg)=(OpTotSq/OpEn)-(AnnualData(XReg)**2)
        if (AnnualStDev(XReg).GT.0) then
          AnnualStDev(XReg)=sqrt(AnnualStDev(XReg))
        else
          AnnualStDev(XReg)=0
        end if
      end if
    end if
  end do
end do

print*, "  > Dumping to twin .agg files..." 

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

GenericFile  = CommFile(XComm,XExec)		! generate region filepaths with correct directory
SuffixStart  = index(GenericFile,".agg")
GenericFile  = GenericFile(1:(SuffixStart-1))
SaveSuffix   = GetFileSuffix (GenericFile)
SuffixStart  = index(GenericFile,SaveSuffix)
call CheckVariSuffix (SaveSuffix,SaveVariName,FileFactor)

GenericFile=CommFile(XComm,XExec) ; LineFormat='(a20,12f8.1,5f8.1)'
GivenFile = GenericFile(1:(SuffixStart-1)) // ".av" // GenericFile(SuffixStart:80)
call SaveAgg (GivenFile,RegNames,Monthly=MonthData,Seasonal=SeasonData,Annual=AnnualData,&
			CallLineFormat=LineFormat,NoResponse=1,ExtraHeads=1, &
			CallVariName=SaveVariName,CallTimeName="mean                ")

GivenFile = GenericFile(1:(SuffixStart-1)) // ".sd" // GenericFile(SuffixStart:80)
LineFormat='(a20,12f8.1,5f8.1)'
call SaveAgg (GivenFile,RegNames,Monthly=MonthStDev,Seasonal=SeasonStDev,Annual=AnnualStDev,&
			CallLineFormat=LineFormat,NoResponse=1,ExtraHeads=1, &
			CallVariName=SaveVariName,CallTimeName="stdev               ")

deallocate (MonthData,SeasonData,AnnualData,RegNames, &
	    MonthStDev,SeasonStDev,AnnualStDev, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: DumpAgg: Deallocation failure: v2 #####"

end subroutine DumpAgg

!*******************************************************************************
! command = 49

subroutine CountOccurrences

call PointOperA			! before we load up into regions, we point OperA towards the correct grim

call RegSelect (nint(MissVal), ExeN, WyeN, (ExeN*WyeN), MapIDLReg, RegSizes, RegNames, &
			RegTitle, RegN, CommInfo(XComm,1))
call GetWeightArrays

allocate (GloSlice(RegN),   stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CountOccurrences: Allocation failure: Gloslice #####"
GloSlice = 0.0		! careful - do not inintialise to MissVal

do XExe = 1, ExeN					! iterate around grid
  do XWye = 1, WyeN					! if the grid cell belongs to a region
    if (MapIDLReg(XExe,XWye).NE.MissVal) then		! then check if the grid cell is within the grim domain
      if (Grid(XExe,XWye).NE.MissVal.AND.RegWeights(MapIDLReg(XExe,XWye)).NE.MissVal) then		
       if (GloSlice(MapIDLReg(XExe,XWye)).NE.MissVal) then
        do XYear = 1, YearN
         do XMonth = 1, MonthN
            if (OperA(XYear,XMonth,Grid(XExe,XWye)).EQ.MissVal.OR. &
            	OperA(XYear,XMonth,Grid(XExe,XWye)).EQ.CommKay(XComm)) &	
                GloSlice(MapIDLReg(XExe,XWye)) = GloSlice(MapIDLReg(XExe,XWye)) + GridWeights(XExe,XWye)
         end do
        end do
       end if
      else						
       GloSlice(MapIDLReg(XExe,XWye)) = MissVal
      end if
    end if
  end do
end do

do XReg = 1, RegN
  if (GloSlice(XReg).NE.MissVal.AND.RegWeights(XReg).NE.MissVal.AND.RegWeights(XReg).GT.0) then
  	GloSlice(XReg) = GloSlice(XReg) / (RegWeights(XReg)*12.0)
  else
        GloSlice(XReg) = MissVal
  end if
end do

nullify (OperA)

deallocate (MapIDLReg, RegSizes, GridWeights, RegWeights, RegBoreal, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CountOccurrences: Deallocation failure: v1 #####"

call SaveGlo (ExeN,WyeN,BoxN,GridFilePath,CommFile(XComm,XExec),CommInfo(XComm,XExec),GloSlice,Grid)

deallocate (GloSlice,RegNames,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CountOccurrences: Deallocation failure: v2 #####"

end subroutine CountOccurrences

!*******************************************************************************
! command = 51

subroutine InterpolateMissing

call PointOperA			
call PointOperB

OperA = OperB

TotA=0 ; TotB=0
do XYear = 1, YearN
 do XMonth = 1, MonthN
  do XExe = 1, ExeN
    do XWye = 1, WyeN
      if (Grid(XExe,XWye).NE.MissVal) then
        if (OperB(XYear,XMonth,Grid(XExe,XWye)).EQ.MissVal) then
          OpTot = 0.0 ; OpEn = 0.0

          do XExeNear = XExe-1,XExe+1
            do XWyeNear = XWye-1,XWye+1
              if (XExeNear.GE.1.AND.XExeNear.LE.ExeN.AND.XWyeNear.GE.1.AND.XWyeNear.LE.WyeN) then
                if (Grid(XExeNear,XWyeNear).NE.MissVal) then
                  if (OperB(XYear,XMonth,Grid(XExeNear,XWyeNear)).NE.MissVal) then
                    OpTot = OpTot + OperB(XYear,XMonth,Grid(XExeNear,XWyeNear))
                    OpEn  = OpEn  + 1.0
                  end if
                end if
              end if
            end do
          end do
          
          if (OpEn.GE.1) then
            OperA(XYear,XMonth,Grid(XExe,XWye)) = OpTot / OpEn
            TotA = TotA + 1
          else
            TotB = TotB + 1
          end if
        end if
      end if
    end do
  end do
 end do
end do

nullify (OperA,OperB)			

end subroutine InterpolateMissing

!*******************************************************************************
! command = 52

subroutine MakeKoeppen

write (99,*), "@@@ point to arrays"		! @@@@@@@@@@@@@@@@@@@@@@@@
call PointOperA			
call PointOperB
write (99,*), "@@@ call RegSelect"		! @@@@@@@@@@@@@@@@@@@@@@@@
call RegSelect (nint(MissVal), ExeN, WyeN, (ExeN*WyeN), MapIDLReg, RegSizes, RegNames, &
			RegTitle, RegN, CommInfo(XComm,1))
write (99,*), "@@@ call GetWeightArrays"	! @@@@@@@@@@@@@@@@@@@@@@@@
call GetWeightArrays
write (99,*), "@@@ make allocations"		! @@@@@@@@@@@@@@@@@@@@@@@@

allocate (MonthTemp(RegN,MonthN),MonthPrec(RegN,MonthN), &
	  GloSlice(RegN),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: MakeKoeppen: Allocation failure."
MonthTemp=0.0 ; MonthPrec=0.0 ; GloSlice=MissVal ; TotA=0 ; TotB=0

XYear = CommOp(XComm,5) - YearAD(1) + 1
  
write (99,*), "@@@ fill MonthData with area means"	! @@@@@@@@@@@@@@@@@@@@@@@@
do XMonth=1,MonthN
 do XExe = 1, ExeN		! fill MonthData with totals			
  do XWye = 1, WyeN					! iterate around grid
    XReg=MapIDLReg(XExe,XWye)		! if the grid cell belongs to a region
    if (XReg.NE.MissVal) then		! then check if the grid cell is within the grim domain
      if (Grid(XExe,XWye).NE.MissVal.AND.RegWeights(XReg).NE.MissVal) then		
        if (OperA(XYear,XMonth,Grid(XExe,XWye)).NE.MissVal.AND. &
	    OperB(XYear,XMonth,Grid(XExe,XWye)).NE.MissVal) then	! check T/P = non-missing
          						! check region = OK thus far
	  if (MonthTemp(XReg,XMonth).NE.MissVal.AND.MonthPrec(XReg,XMonth).NE.MissVal) then
	    MonthTemp(XReg,XMonth) = MonthTemp(XReg,XMonth) + &
		(OperA(XYear,XMonth,Grid(XExe,XWye)) * GridWeights(XExe,XWye))
	    MonthPrec(XReg,XMonth) = MonthPrec(XReg,XMonth) + &
		(OperB(XYear,XMonth,Grid(XExe,XWye)) * GridWeights(XExe,XWye))
          end if
	else						! datum=missing, so make region missing
          MonthTemp(XReg,XMonth)=MissVal ; MonthPrec(XReg,XMonth)=MissVal
        end if
      else						! cell=missing, so make region missing			
        MonthTemp(XReg,XMonth)=MissVal ; MonthPrec(XReg,XMonth)=MissVal
      end if
    end if
  end do
 end do

 do XReg=1,RegN			! convert area totals to area means
   if (RegWeights(XReg).GT.0.0) then
     if (MonthTemp(XReg,XMonth).NE.MissVal.AND.MonthPrec(XReg,XMonth).NE.MissVal) then
          MonthTemp(XReg,XMonth) = MonthTemp(XReg,XMonth) / RegWeights(XReg)
          MonthPrec(XReg,XMonth) = MonthPrec(XReg,XMonth) / RegWeights(XReg)
     end if
   else
     MonthTemp(XReg,XMonth)=MissVal ; MonthPrec(XReg,XMonth)=MissVal
   end if
 end do
end do

write (99,*), "@@@ calculate Koeppen class"		! @@@@@@@@@@@@@@@@@@@@@@@@
do XReg=1,RegN			! calculate Koeppen class for each region
  do XMonth=1,MonthN
    BoxTemp(XMonth)=MonthTemp(XReg,XMonth) ; BoxPrec(XMonth)=MonthPrec(XReg,XMonth)
  end do
  Boreal=.TRUE. ; if (RegBoreal(XReg).LT.0.0) Boreal=.FALSE.
  
  call GetKoeppen (BoxTemp,BoxPrec,Boreal,KopNumb,KopAcro,KopName)
  
  if (KopNumb.NE.MissVal) then
    TotA=TotA+1 
  else 
    TotB=TotB+1
  end if
  
  GloSlice(XReg)=KopNumb
end do

InfoLine="Koeppen classes"
write (99,*), "@@@ call SaveGlo"		! @@@@@@@@@@@@@@@@@@@@@@@@
call SaveGlo (ExeN,WyeN,BoxN,GridFilePath,CommFile(XComm,XExec),InfoLine,GloSlice,Grid)
write (99,*), "@@@ done SaveGlo"		! @@@@@@@@@@@@@@@@@@@@@@@@

deallocate (MonthTemp,MonthPrec,GloSlice,MapIDLReg,RegSizes,RegNames, &
		GridWeights,RegWeights,RegBoreal,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: MakeKoeppen: Deallocation failure."

nullify (OperA,OperB)			
write (99,*), "@@@ done MakeKoeppen"		! @@@@@@@@@@@@@@@@@@@@@@@@

end subroutine MakeKoeppen

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

subroutine GetWeightArrays

if (CommOp(XComm,4).EQ.2) then
  call GetGloWeights (ExeN,WyeN,GridWeights)
else
  allocate (GridWeights(ExeN,WyeN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: GetWeightArrays: Allocation failure: GridWeights #####"
  GridWeights = 1.0
end if

allocate (RegWeights(RegN), RegBoreal(RegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetGridWeights: Allocation failure: RegWeights #####"
RegWeights = 0.0 ; RegBoreal = 0.0	! careful - do not initialise to MissVal

BoxLat=(Bounds(4)-Bounds(3))/WyeN

do XWye = 1, WyeN
  NorthSouth=1.0
  if (((real(XWye)*BoxLat)+Bounds(3)).LT.0) NorthSouth=-1.0
  
  do XExe = 1, ExeN
    XReg=MapIDLReg(XExe,XWye)
    if (XReg.NE.MissVal.AND.Grid(XExe,XWye).NE.MissVal) then
      if (RegWeights(XReg).NE.MissVal.AND.GridWeights(XExe,XWye).NE.MissVal) then
        RegWeights(XReg) = RegWeights(XReg) + GridWeights(XExe,XWye)
        RegBoreal (XReg) = RegBoreal(XReg) + (GridWeights(XExe,XWye) * NorthSouth)
      else
        RegWeights(XReg) = MissVal ; RegBoreal(XReg) = MissVal
      end if
    end if
  end do
end do

end subroutine GetWeightArrays

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

subroutine PointOperA

if (CommOp(XComm,2).EQ.1) OperA => Grim1
if (CommOp(XComm,2).EQ.2) OperA => Grim2
if (CommOp(XComm,2).EQ.3) OperA => Grim3
if (CommOp(XComm,2).EQ.4) OperA => Grim4
if (CommOp(XComm,2).EQ.5) OperA => Grim5
if (CommOp(XComm,2).EQ.6) OperA => Grim6

end subroutine PointOperA

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

subroutine PointOperB

if (CommOp(XComm,3).EQ.1) OperB => Grim1
if (CommOp(XComm,3).EQ.2) OperB => Grim2
if (CommOp(XComm,3).EQ.3) OperB => Grim3
if (CommOp(XComm,3).EQ.4) OperB => Grim4
if (CommOp(XComm,3).EQ.5) OperB => Grim5
if (CommOp(XComm,3).EQ.6) OperB => Grim6

end subroutine PointOperB

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

subroutine PointOperC

if (CommOp(XComm,5).EQ.1) OperC => Grim1
if (CommOp(XComm,5).EQ.2) OperC => Grim2
if (CommOp(XComm,5).EQ.3) OperC => Grim3
if (CommOp(XComm,5).EQ.4) OperC => Grim4
if (CommOp(XComm,5).EQ.5) OperC => Grim5
if (CommOp(XComm,5).EQ.6) OperC => Grim6

end subroutine PointOperC

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

subroutine PointOperD

if (CommOp(XComm,4).EQ.1) OperD => Grim1
if (CommOp(XComm,4).EQ.2) OperD => Grim2
if (CommOp(XComm,4).EQ.3) OperD => Grim3
if (CommOp(XComm,4).EQ.4) OperD => Grim4
if (CommOp(XComm,4).EQ.5) OperD => Grim5
if (CommOp(XComm,4).EQ.6) OperD => Grim6

end subroutine PointOperD

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

subroutine PointOperK

if (CommKay(XComm).EQ.1) OperK => Grim1
if (CommKay(XComm).EQ.2) OperK => Grim2
if (CommKay(XComm).EQ.3) OperK => Grim3
if (CommKay(XComm).EQ.4) OperK => Grim4
if (CommKay(XComm).EQ.5) OperK => Grim5
if (CommKay(XComm).EQ.6) OperK => Grim6

end subroutine PointOperK

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

subroutine GridMismatch

if (size(FileGrid,1).NE.size(Grid,1)) print*, "  > ##### ERROR: loaded Grid mismatch: X"
if (size(FileGrid,2).NE.size(Grid,2)) print*, "  > ##### ERROR: loaded Grid mismatch: Y"

CheckGrid = 0
do XExe = 1, ExeN
  do XWye = 1, WyeN
    if (FileGrid(XExe,XWye).NE.Grid(XExe,XWye)) then
      CheckGrid = CheckGrid + 1
      write (99,"(4i10)"), XExe,XWye,Grid(XExe,XWye),FileGrid(XExe,XWye)
    end if
  end do
end do
if (CheckGrid.NE.0) print*, "  > ##### ERROR: Grid mismatches: ", CheckGrid

end subroutine GridMismatch

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

subroutine Finalise

deallocate (CommOp,CommKay,CommFile,CommInfo, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Finalise: Deallocation failure: Comm* #####"

deallocate (Grid,YearAD,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Finalise: Deallocation failure: Basic #####"

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

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

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

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

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

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

print*

close (99)

end subroutine Finalise

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

end program SeqGrim