! setupseqgrim.f90
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
! 	-o ./../grim/setupseqgrim time.f90 filenames.f90 grimfiles.f90 
!	glofiles.f90 initialmod.f90 grid.f90 ./../grim/setupseqgrim.f90 
!	2> /tyn1/tim/scratch/stderr.txt
! written by Tim Mitchell on 02.04.01
! last modified on 23.10.01
! this program constructs the specification file for use in SeqGrim

program SetupSeqGrim

use Time
use FileNames
use GrimFiles
use GloFiles
use InitialMod
use Grid

implicit none

real, pointer, dimension (:)			:: CommKay

integer, pointer, dimension (:,:)		:: CommOp, RefGrid, RSMapIDLReg
integer, pointer, dimension (:) 		:: RSRegSizes	

character (len=80), pointer, dimension (:,:)	:: CommFile, CommInfo
character (len=20), pointer, dimension (:) 	:: RSRegNames

real, dimension (4)				:: RefBounds

real, parameter 		:: MissVal = -999.0

real :: ThisLon,ThisLat

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

integer :: RSGridChosen,RSGridLongN,RSGridLatN,RSGridDataN			
integer :: AllocStat, ReadStatus
integer :: YearN, MonthN, BoxN, ExeN, WyeN, ExecN, ArrayN, CommN, RSRegN
integer :: XYear, XMonth, XBox, XExe, XWye, XExec, XArray, XComm, XElement, XFileLoad, XBound
integer :: SubLen, FileSubBeg, InfoSubBeg, OrigSubBeg, Store
integer :: RefBoxN, YearAD0, YearAD1
integer :: QHowGB, QWeight

character (len=80) :: SpecFile, FileInfo, GivenFile, CheckFile, GivenInfo, OrigSub, SpecSub, TextString
character (len=80) :: RSGridFile			
character (len=40) :: RSRegTitle			
character (len=10) :: RSGridTitle			
character (len= 4) :: FileSuffix, CheckSuffix, SaveSuffix

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

open  (99,file="./../../../scratch/log-susqg.dat",status="replace",action="write")

call Initialise
call FirstExec
if (ExecN.GT.1) call AutoSpec
call SaveSpecFile
call Finalise

close (99)

contains

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

subroutine Initialise

print*
print*, "  > ***** SetUpSeqGrim: prepares spec file for SeqGrim *****"
print*

print*, "  > The array structure is common throughout."
print*, "  > Submit a grim that contains the appropriate grid."
call GrabGrid (1,RefGrid,RefBounds,RefBoxN,Quiet=1)

ExeN = size (RefGrid,1) ; WyeN = size (RefGrid,2) ; BoxN = RefBoxN

print*, "  > Enter the first, last years AD: "
do
	read (*,*,iostat=ReadStatus), YearAD0, YearAD1
	if (ReadStatus.LE.0.AND.YearAD1.GE.YearAD0) exit
end do
YearN  = YearAD1 - YearAD0 + 1
MonthN = 12

print*, "  > Enter the number of executions of the commands: "
do
	read (*,*,iostat=ReadStatus), ExecN
	if (ReadStatus.LE.0.AND.ExecN.GE.1) exit
end do

print*, "  > Enter the number of commands per execution: "
do
	read (*,*,iostat=ReadStatus), CommN
	if (ReadStatus.LE.0.AND.CommN.GE.1) exit
end do

print*, "  > Enter the number of arrays required (1-6): "
do
	read (*,*,iostat=ReadStatus), ArrayN
	if (ReadStatus.LE.0.AND.ArrayN.GE.1.AND.ArrayN.LE.6) exit
end do

print*, "  > Enter the filepath of the operations file (.ops): "
print*, "  >   (Place in ./../../../code/f90/grim/_ref/???.ops)"
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SpecFile = SavePath (GivenFile,".ops")
! SpecFile = trim(SpecFile) // ".X"

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

end subroutine Initialise

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

subroutine FirstExec

print*
if (ExecN.EQ.1) print*, "  > We iterate through the execution, specifying each command. "
if (ExecN.GT.1) print*, "  > We iterate through the first execution, specifying each command. "
print*, "  > Load: -1:A->reset,0:->A,1:yr->A,2:lin->A,3:ann->A,4:->A(per)"
print*, "  > Aver: 10:A=B(perav),11:A=B(gloav),12:A=B(annav),13:A=A(anom),14:A=Gauss(B)"
print*, "  > Also: 15:A=B(c=>d),16:where(A=C)B=k,17:where(A[rel=e]C)B=D,18:where(A=k)B=C"
print*, "  >       19:where(A[rel=c]k)B=k,20:LSR(y=ax+b)D=AC+B,21:A=monlen,22:A=B(trunc)"
print*, "  > Oper: 30:A=B.k,31:A=B.C,32:A=k,33=A.external-year"
print*, "  > Save: 43:A->grim,44:A->period,45:A->.glo,"
print*, "  >       46:A->.per,47:A->reg->.per,48:A->.agg,49:count(A=k)->.glo"
print*, "  > Comm: 50:mean/sum,51:interpol,52:Koeppen"

do XComm = 1, CommN
  print "(a,i3)", "   > Command (-10=list): ", XComm  
  do
	read (*,*,iostat=ReadStatus), CommOp(XComm,1)
        if (CommOp(XComm,1).EQ.-10) then
	 print*, "  > Load: -1:A->reset,0:->A,1:yr->A,2:lin->A,3:ann->A,4:->A(per)"
	 print*, "  > Aver: 10:A=B(perav),11:A=B(gloav),12:A=B(annav),13:A=A(anom),14:A=Gauss(B)"
	 print*, "  > Also: 15:A=B(c=>d),16:where(A=C)B=k,17:where(A[rel=e]C)B=D,18:where(A=k)B=C,"
	 print*, "  >       19:where(A[rel=c]k)B=k,20:LSR(y=ax+b)D=AC+B,21:A=monlen,22:A=B(trunc)"
	 print*, "  > Oper: 30:A=B.k,31:A=B.C,32:A=k,33=A.external-year"
	 print*, "  > Save: 43:A->grim,44:A->period,45:A->.glo,"
	 print*, "  >       46:A->.per,47:A->reg->.per,48:A->.agg,49:count(A=k)->.glo"
	 print*, "  > Comm: 50:mean/sum,51:interpol,52:Koeppen"
	end if
	if (ReadStatus.LE.0.AND.CommOp(XComm,1).GE.-1.AND.CommOp(XComm,1).LE.52) exit
  end do
  
  if      (CommOp(XComm,1).EQ.-1) then			! reset A
    call SelectA
  else if (CommOp(XComm,1).EQ. 0) then			! load A - any period overlapping with main array
    call SelectA					! any non-overlapping period is left unchanged
!    call SelectFL					! to initialise main array first, use option (-1)
    if (ExecN.EQ.1) print*, "  > Enter the filepath: " 						
    if (ExecN.GT.1) print*, "  > Enter the filepath (including substring to vary later): " 	
    call SelectFP    
    call ReviewCall (GivenFile,"    ",CheckFile,CheckSuffix,1)
    CommFile (XComm,1) = CheckFile
  else if (CommOp(XComm,1).EQ. 1) then			! load A single year
    call SelectA
    if (ExecN.EQ.1) print*, "  > Enter the filepath: " 						
    if (ExecN.GT.1) print*, "  > Enter the filepath (including substring to vary later): " 	
    call SelectFP
    call ReviewCall (GivenFile,"    ",CheckFile,CheckSuffix,1)
    CommFile (XComm,1) = CheckFile
  else if (CommOp(XComm,1).EQ. 2) then			! load A .lin
    call SelectA
    call SelectPP
    if (ExecN.EQ.1) print*, "  > Enter the filepath: " 						
    if (ExecN.GT.1) print*, "  > Enter the filepath (including substring to vary later): " 	
    call SelectFP    
    CommFile (XComm,1) = LoadPath (GivenFile,".lin")
  else if (CommOp(XComm,1).EQ. 3) then			! load A .ann (ID column)
    call SelectA
    call SelectCL
    if (ExecN.EQ.1) print*, "  > Enter the filepath: " 						
    if (ExecN.GT.1) print*, "  > Enter the filepath (including substring to vary later): " 	
    call SelectFP    
    CommFile (XComm,1) = LoadPath (GivenFile,".ann")
  else if (CommOp(XComm,1).EQ. 4) then			! load grim into A (specific period)
    call SelectA
    print*, "  > Load grim file into a specific period."					
    call SelectFL					
    if (ExecN.EQ.1) print*, "  > Enter the filepath: " 						
    if (ExecN.GT.1) print*, "  > Enter the filepath (including substring to vary later): " 	
    call SelectFP    
    call ReviewCall (GivenFile,"    ",CheckFile,CheckSuffix,1)
    CommFile (XComm,1) = CheckFile
  else if (CommOp(XComm,1).EQ.10) then			! A = av(B)
    call SelectA
    call SelectB
    call SelectFL
  else if (CommOp(XComm,1).EQ.11) then			! glo av
    call SelectA
    call SelectB
  else if (CommOp(XComm,1).EQ.12) then			! ann av
    call SelectA
    call SelectB
  else if (CommOp(XComm,1).EQ.13) then			! A = anom(A)
    call SelectA
    call SelectFL
    call SelectRM
  else if (CommOp(XComm,1).EQ.14) then			! A = Gauss(B)
    call SelectA
    call SelectB
    call SelectPQ
  else if (CommOp(XComm,1).EQ.15) then			! A = (B) change c => d
    call SelectA
    call SelectB
    call SelectCD
  else if (CommOp(XComm,1).EQ.16) then			! where A=c, B=k
    call SelectA
    call SelectB
    call SelectCK
  else if (CommOp(XComm,1).EQ.17) then			! where A[rel=e]C, B=D
    call SelectA
    call SelectB
    call SelectC
    call SelectD
    call SelectRL
  else if (CommOp(XComm,1).EQ.18) then			! where A=k, B=C
    call SelectA
    call SelectB
    call SelectCK 	! fills 'D' (col4) and k
  else if (CommOp(XComm,1).EQ.19) then			! where A[rel=c]k, B=k
    call SelectA
    call SelectB
    call SelectIG
    call SelectRL
    CommOp(XComm,5) = nint(CommKay(XComm))
    call SelectK
  else if (CommOp(XComm,1).EQ.20) then			! LSR - look at SelectB before changing '20'
    call SelectA
    print*, "  > If b=-999 use y=ax (D=AC), else use y=ax+b (D=AC+B)"
    call SelectB
    print*, "  > If 'y' and 'x' are the same, we regress 'y' against time."
    call SelectC
    call SelectD
    if (CommOp(XComm,3).NE.MissVal) call SelectKCorr
  else if (CommOp(XComm,1).EQ.21) then			! month lengths
    call SelectA
  else if (CommOp(XComm,1).EQ.22) then			! truncate
    call SelectA
    call SelectB
    call SelectTL
  else if (CommOp(XComm,1).EQ.30) then			! A = B.k
    call SelectA
    call SelectB
    call SelectOK
  else if (CommOp(XComm,1).EQ.31) then			! A = B.C
    call SelectA
    call SelectB
    call SelectOP
    call SelectC
  else if (CommOp(XComm,1).EQ.32) then			! A = k
    call SelectA
    call SelectK
  else if (CommOp(XComm,1).EQ.33) then			! A = B.external-single-year-grim
    call SelectA
    call SelectB
    call SelectOP
    if (ExecN.EQ.1) print*, "  > Enter the filepath: " 						
    if (ExecN.GT.1) print*, "  > Enter the filepath (including substring to vary later): " 	
    call SelectFP
    call ReviewCall (GivenFile,"    ",CheckFile,CheckSuffix,1)
    CommFile (XComm,1) = CheckFile
  else if (CommOp(XComm,1).EQ.43) then			! save A
    call SelectA
    if (ExecN.EQ.1) print*, "  > Enter the filepath: " 						
    if (ExecN.GT.1) print*, "  > Enter the filepath (including substring to vary later): " 	
    call SelectFP
    call ReviewCall (GivenFile,"    ",CheckFile,CheckSuffix,2)
    CommFile (XComm,1) = CheckFile
    print*, "  > Enter the information line: "
    call SelectIL
  else if (CommOp(XComm,1).EQ.44) then			! save A period
    call SelectA
    call SelectFL
    if (ExecN.EQ.1) print*, "  > Enter the filepath: " 						
    if (ExecN.GT.1) print*, "  > Enter the filepath (including substring to vary later): " 	
    call SelectFP
    call ReviewCall (GivenFile,"    ",CheckFile,CheckSuffix,2)
    CommFile (XComm,1) = CheckFile
    print*, "  > Enter the information line: "
    call SelectIL
  else if (CommOp(XComm,1).EQ.45) then			! save A .glo
    call SelectA
    call SelectOY
    call SelectAV
    TextString = "  > Enter the .glo file"
    if (ExecN.GT.1) TextString = trim(TextString) // " (including substring to vary later)" 	
    if (CommOp(XComm,4).EQ.2) TextString = trim(TextString) // " (including 'MAM')"
    if (CommOp(XComm,4).EQ.3) TextString = trim(TextString) // " (including '01')"
    TextString = trim(TextString) // ":"
    print*, trim(TextString)
    call SelectFP    
    CommFile (XComm,1) = SavePath (GivenFile,".glo")
    if (CommOp(XComm,4).EQ.1) print*, "  > Enter the information line: "
    if (CommOp(XComm,4).EQ.2) print*, "  > Enter the information line (including 'MAM'): "
    if (CommOp(XComm,4).EQ.3) print*, "  > Enter the information line (including '01'): "
    call SelectIL
  else if (CommOp(XComm,1).EQ.46) then			! save A .per
    call SelectA
    call SelectGB
    print*, "  > The save file must have a variable suffix: we add .per automatically."
    if (ExecN.EQ.1) print*, "  > Enter the file: " 						
    if (ExecN.GT.1) print*, "  > Enter the file (with substring to vary): " 	
    do
      call SelectFP
      SaveSuffix = GetFileSuffix (GivenFile)
      if (SaveSuffix.EQ."    ") print*, "  > Try again."
      if (SaveSuffix.NE."    ") exit
    end do
    GivenFile = trim(GivenFile) // ".per"
    CommFile (XComm,1) = SavePath (GivenFile,".per")
  else if (CommOp(XComm,1).EQ.47) then			! A --> region set --> save each region to .per
    call SelectA
    print*, "  > Identify the grid and region set for which to calc regional values."
    call SelectRS
    print*, "  > Weight by box=1 (=1), or by area (=2) ?"
    call SelectWE
    print*, "  > Give the save file (maxlen=56) a variable suffix: reg-name + .per = auto."
    if (ExecN.EQ.1) print*, "  > Enter the file: " 						
    if (ExecN.GT.1) print*, "  > Enter the file (including substring to vary later): " 	
    do
      call SelectFP
      SaveSuffix = GetFileSuffix (GivenFile)
      if (SaveSuffix.EQ."    ") print*, "  > Try again."
      if (len_trim(GivenFile).GT.56) print*, "  >  Max length of filepath is 56 characters. Try again."
      if (len_trim(GivenFile).LE.56.AND.SaveSuffix.NE."    ") exit
    end do
    GivenFile = trim(GivenFile) // ".per"
    CommFile (XComm,1) = SavePath (GivenFile,".per")
    print*, "  > WARNING: You may need to alter the optional arguments in :"
    print*, "  >          seqgrim.f90 --> DumpDotGlo --> SavePer"
  else if (CommOp(XComm,1).EQ.48) then			! save A .agg (both mean and stdev)
    call SelectA
    call SelectFL
    CommOp(XComm,3)=CommOp(XComm,4)	! change period to 3,5 not 4,5
    print*, "  > Identify the grid and region set for which to calc regional values."
    call SelectRS
    print*, "  > Weight by box=1 (=1), or by area (=2) ?"
    call SelectWE			! weight goes in 4
    print*, "  > Give the save file a variable suffix: .agg = auto."
    if (ExecN.EQ.1) print*, "  > Enter the file: " 						
    if (ExecN.GT.1) print*, "  > Enter the file (including substring to vary later): " 	
    do
      call SelectFP
      SaveSuffix = GetFileSuffix (GivenFile)
      if (SaveSuffix.EQ."    ") print*, "  > Try again."
      if (SaveSuffix.NE."    ") exit
    end do
    GivenFile = trim(GivenFile) // ".agg"
    CommFile (XComm,1) = SavePath (GivenFile,".agg")
    print*, "  > WARNING: You may need to alter the optional arguments in :"
    print*, "  >          seqgrim.f90 --> DumpAgg --> SaveAgg"
  else if (CommOp(XComm,1).EQ.49) then			! count A=k -> .glo
    call SelectA
    call SelectK
    print*, "  > Identify the grid and region set for which to calc regional values."
    call SelectRS
    print*, "  > Weight by box=1 (=1), or by area (=2) ?"
    call SelectWE
    print*, "  > Give the save file a variable suffix: .glo = auto."
    if (ExecN.EQ.1) print*, "  > Enter the file: " 						
    if (ExecN.GT.1) print*, "  > Enter the file (including substring to vary later): " 	
    do
      call SelectFP
      SaveSuffix = GetFileSuffix (GivenFile)
      if (SaveSuffix.EQ."    ") print*, "  > Try again."
      if (SaveSuffix.NE."    ") exit
    end do
    GivenFile = trim(GivenFile) // ".glo"
    CommFile (XComm,1) = SavePath (GivenFile,".glo")
  else if (CommOp(XComm,1).EQ.50) then			! override the auto mean/sum feature
    print*, "  > Override the auto. Select mean (=0) or sum (=2)."
    do
	  read (*,*,iostat=ReadStatus), CommOp(XComm,2)
	  if (ReadStatus.LE.0.AND.CommOp(XComm,2).GE.0.AND.CommOp(XComm,2).LE.2) exit
    end do
  else if (CommOp(XComm,1).EQ.51) then			! spatially interpolate missing A=B
    call SelectA
    call SelectB
  else if (CommOp(XComm,1).EQ.52) then			! calc Koeppen for regions -> .glo
    print*, "  > A=temp(oC) B=prec(mm)"
    call SelectA
    call SelectB  ; Store=CommOp(XComm,3)
    call SelectOY ; CommOp(XComm,5)=CommOp(XComm,3) ; CommOp(XComm,3)=Store
    print*, "  > Identify the grid and region set for which to calc regional values."
    call SelectRS
    print*, "  > Weight by box=1 (=1), or by area (=2) ?"
    call SelectWE
    print*, "  > Enter the Koeppen .glo file to save." 
    call SelectFP    
    CommFile (XComm,1) = SavePath (GivenFile,".glo")
  end if
  
  print*
end do

end subroutine FirstExec

!*******************************************************************************
! select A

subroutine SelectA

if (ArrayN.GT.1) then
  print "(a,i1,a2)", "   > Select A  (1...", ArrayN, "):"	
  do
	read (*,*,iostat=ReadStatus), CommOp(XComm,2)
	if (ReadStatus.LE.0.AND.CommOp(XComm,2).GE.1.AND.CommOp(XComm,2).LE.ArrayN) exit
  end do
else
  CommOp(XComm,2) = 1
end if

end subroutine SelectA

!*******************************************************************************
! select averaging period --> 4

subroutine SelectAV

print*, "  > Select the averaging period (1=annual,2=seasonal,3=monthly):"
do
	  read (*,*,iostat=ReadStatus), CommOp(XComm,4)
	  if (ReadStatus.LE.0.AND.CommOp(XComm,4).GE.1.AND.CommOp(XComm,4).LE.3) exit
end do
    
end subroutine SelectAV

!*******************************************************************************
! select B

subroutine SelectB

if (ArrayN.GT.1) then
  print "(a,i1,a2)", "   > Select B  (1...", ArrayN, "):"	
  do
	read (*,*,iostat=ReadStatus), CommOp(XComm,3)
	if (ReadStatus.LE.0.AND.CommOp(XComm,1).EQ.20.AND. CommOp(XComm,3)       .EQ.MissVal) exit
	if (ReadStatus.LE.0.AND.CommOp(XComm,3).GE.1 .AND. CommOp(XComm,3)       .LE.ArrayN)  exit
  end do
else
  CommOp(XComm,3) = 1
end if

end subroutine SelectB

!*******************************************************************************
! select C

subroutine SelectC

if (ArrayN.GT.1) then
  print "(a,i1,a2)", "   > Select C  (1...", ArrayN, "):"	
  do
	read (*,*,iostat=ReadStatus), CommOp(XComm,5)	
	if (ReadStatus.LE.0.AND.CommOp(XComm,5).GE.1.AND.CommOp(XComm,5).LE.ArrayN) exit
  end do
else
  CommOp(XComm,5) = 1
end if

end subroutine SelectC

!*******************************************************************************
! select C and D (i.e. constants, but integers)

subroutine SelectCD

print*, "  > Select 'c' and 'd':"	
do
	read (*,*,iostat=ReadStatus), CommOp(XComm,4),CommOp(XComm,5)	
	if (ReadStatus.LE.0) exit
end do

end subroutine SelectCD

!*******************************************************************************
! select C and k (i.e. constants, but C=integer)

subroutine SelectCK

print*, "  > Select 'c' and 'k':"	
do
	read (*,*,iostat=ReadStatus), CommOp(XComm,4),CommKay(XComm)	
	if (ReadStatus.LE.0) exit
end do

end subroutine SelectCK

!*******************************************************************************
! select column

subroutine SelectCL

print "(a,i1,a2)", "   > Select the .ann column to load:"	
do
	read (*,*,iostat=ReadStatus), CommOp(XComm,4)	
	if (ReadStatus.LE.0.AND.CommOp(XComm,4).GE.1) exit
end do

end subroutine SelectCL

!*******************************************************************************
! select D

subroutine SelectD

if (ArrayN.GT.1) then
  print "(a,i1,a2)", "   > Select D  (1...", ArrayN, "):"	
  do
	read (*,*,iostat=ReadStatus), CommOp(XComm,4)
	if (ReadStatus.LE.0.AND.CommOp(XComm,4).GE.1.AND.CommOp(XComm,4).LE.ArrayN) exit
  end do
else
  CommOp(XComm,4) = 1
end if

end subroutine SelectD

!*******************************************************************************
! select period --> 4,5

subroutine SelectFL

print*, "  > Select the first, last years AD:"
do
  read (*,*,iostat=ReadStatus), CommOp(XComm,4), CommOp(XComm,5)
  if (CommOp(XComm,4).GT.CommOp(XComm,5).OR.CommOp(XComm,4).LT.YearAD0.OR.CommOp(XComm,5).GT.YearAD1) then
	    print*, "  > Unacceptable period. Try again."
	    ReadStatus = 1
  end if 	  
  if (ReadStatus.LE.0) exit
end do
    
end subroutine SelectFL

!*******************************************************************************
! select filepath

subroutine SelectFP

do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
    
end subroutine SelectFP

!*******************************************************************************
! select grid box --> 3

subroutine SelectGB

CommOp(XComm,3) = MissVal

print*, "  > Select grid box: by box no. (=0), grid-ref (=1), or a long/lat (=2) ?"
do
	  read (*,*,iostat=ReadStatus), QHowGB
	  if (ReadStatus.LE.0.AND.QHowGB.GE.0.AND.QHowGB.LE.2) exit
end do

do
  if (QHowGB.EQ.0) then
    print*, "  > Enter the box number:"
    do
	  read (*,*,iostat=ReadStatus), XBox
	  if (ReadStatus.LE.0) exit
    end do
  else if (QHowGB.EQ.1) then    
    print*, "  > Enter the grid-ref (x,y):"
    do
	  read (*,*,iostat=ReadStatus), XExe, XWye
	  if (XExe.GE.1.AND.XExe.LE.ExeN.AND.XWye.GE.1.AND.XWye.LE.WyeN) then
	  	! no prob
	  else
	  	ReadStatus = 1
	  	print*, "  > Outside bounds of grid. Try again."
	  end if
	  if (ReadStatus.LE.0) exit
    end do
  else
    print*, "  > Enter the longitude and latitude: "
    do
	read (*,*,iostat=ReadStatus), ThisLon, ThisLat
	if (ThisLon.LT.RefBounds(1).OR.ThisLon.GT.RefBounds(2) &
			.OR.ThisLat.LT.RefBounds(3).OR.ThisLat.GT.RefBounds(4)) then
	  	ReadStatus = 1
	  	print*, "  > Outside bounds of grid. Try again."
	end if
	if (ReadStatus.LE.0) exit
    end do
    
    XExe = ceiling (real(ExeN) * (ThisLon-RefBounds(1)) / (RefBounds(2)-RefBounds(1)))
    if (XExe.LT.1)    XExe = 1
    if (XExe.GT.ExeN) XExe = ExeN
    
    XWye = ceiling (real(WyeN) * (ThisLat-RefBounds(3)) / (RefBounds(4)-RefBounds(3)))
    if (XWye.LT.1)    XWye = 1
    if (XWye.GT.WyeN) XWye = WyeN
    
    print "(a,2i6)", "   > Grid reference (x,y): ", XExe, XWye    
  end if
  
  if (QHowGB.EQ.0) then
    CommOp(XComm,3) = XBox
  else  
    if (RefGrid(XExe,XWye).EQ.MissVal) print*, "  > This box is set to missing. Try again."
    if (RefGrid(XExe,XWye).NE.MissVal) CommOp(XComm,3) = RefGrid(XExe,XWye)
  end if
  if (CommOp(XComm,3).NE.MissVal) exit
end do
    
end subroutine SelectGB

!*******************************************************************************
! decide whether or not to ignore the missing values

subroutine SelectIG

print*, "  > Observe (=0) or ignore (=1) missing values ?"
do
	read (*,*,iostat=ReadStatus), CommOp(XComm,4)
	if (ReadStatus.LE.0.AND.CommOp(XComm,4).GE.0.AND.CommOp(XComm,4).LE.1) exit
end do      
    
end subroutine SelectIG

!*******************************************************************************
! select info line

subroutine SelectIL

do
	read (*,*,iostat=ReadStatus), CommInfo (XComm,1)
	if (ReadStatus.LE.0.AND.CommInfo (XComm,1).NE."") exit
end do      
    
end subroutine SelectIL

!*******************************************************************************
! select constant

subroutine SelectK

print*, "  > Enter the constant:"
do
	read (*,*,iostat=ReadStatus), CommKay(XComm)
	if (ReadStatus.LE.0) exit
end do      
    
end subroutine SelectK

!*******************************************************************************
! select constant

subroutine SelectKCorr

print*, "  > Decide where (-999=don't) to place the correlation coefficient:"
do
	read (*,*,iostat=ReadStatus), CommKay(XComm)
	if (ReadStatus.LE.0.AND.(CommKay(XComm).EQ.MissVal.OR. &
		(CommKay(XComm).GE.1.AND.CommKay(XComm).LE.ArrayN))) exit
end do      
    
end subroutine SelectKCorr

!*******************************************************************************
! select op and K

subroutine SelectOK

print*, "  > Select op: 1=B/k,2=B*k,3=B-k,4=B+k,"
print*, "               5=B^.5,6=B^k,7=|B|,8=log10(B),9=exp(B),10=ln(B)"	! operation for B.k
do
	read (*,*,iostat=ReadStatus), CommOp(XComm,4)
	if (ReadStatus.LE.0.AND.CommOp(XComm,4).GE.1.AND.CommOp(XComm,4).LE.10) exit
end do
      
if (CommOp(XComm,4).NE.5.AND.CommOp(XComm,4).LT.7) then			! constant
        print*, "  > Enter the constant:"
        do
	  read (*,*,iostat=ReadStatus), CommKay(XComm)
	  if (ReadStatus.LE.0) exit
        end do
end if

end subroutine SelectOK

!*******************************************************************************
! select op only

subroutine SelectOP

print*, "  > Select op (1=B/C,2=B*C,3=B-C,4=B+C):"			! operation for B.C
do
	read (*,*,iostat=ReadStatus), CommOp(XComm,4)
	if (ReadStatus.LE.0.AND.CommOp(XComm,4).GE.1.AND.CommOp(XComm,4).LE.4) exit
end do

end subroutine SelectOP

!*******************************************************************************
! select one year to save --> 3

subroutine SelectOY

if (YearN.GT.1) then
  print*, "  > Select the year AD:"
  do
	  read (*,*,iostat=ReadStatus), CommOp(XComm,3)
	  if (ReadStatus.LE.0.AND.CommOp(XComm,3).GE.YearAD0.AND.CommOp(XComm,3).LE.YearAD1) exit
  end do
else
  CommOp(XComm,3)=YearAD0
end if
    
end subroutine SelectOY

!*******************************************************************************
! select precise period --> 3

subroutine SelectPP

print "(a,i4,a1,i4,a)", "   > Is the file period ", YearAD0, "-", YearAD1, " ? (1=yes,2=no)"
do
	  read (*,*,iostat=ReadStatus), CommOp(XComm,3)
	  if (ReadStatus.LE.0.AND.CommOp(XComm,3).GE.1.AND.CommOp(XComm,3).LE.2) exit
end do
    
end subroutine SelectPP

!*******************************************************************************
! select smoothing period and QRetainMiss

subroutine SelectPQ

print*, "  > Enter the period at which to smooth (years): "
do
	  read (*,*,iostat=ReadStatus), CommOp(XComm,4)
	  if (ReadStatus.LE.0.AND.CommOp(XComm,4).GE.1) exit
end do
    
print*, "  > Replace (=1) or retain (=2) missing values ? "
do
	  read (*,*,iostat=ReadStatus), CommOp(XComm,5)
	  if (ReadStatus.LE.0.AND.CommOp(XComm,5).GE.1.AND.CommOp(XComm,5).LE.2) exit
end do
    
end subroutine SelectPQ

!*******************************************************************************
! select relationship to use

subroutine SelectRL

print*, "  > Enter the relationship (-3=NE,-2=LT,-1=LE,0=EQ,1=GE,2=GT): "
do
	  read (*,*,iostat=ReadStatus), CommKay(XComm)
	  if (ReadStatus.LE.0.AND.CommKay(XComm).GE.-3.AND.CommKay(XComm).LE.2) exit
end do
    
end subroutine SelectRL

!*******************************************************************************
! select whether to (=1, not=0) retain values where mean is missing

subroutine SelectRM

print*, "  > Retain values where mean is missing ? (0=no,1=yes)"
do
	  read (*,*,iostat=ReadStatus), CommOp(XComm,3)
	  if (ReadStatus.LE.0.AND.CommOp(XComm,3).GE.0.AND.CommOp(XComm,3).LE.1) exit
end do
    
end subroutine SelectRM

!*******************************************************************************
! select region set --> CommInfo
! note that if the grim file is on the h2grid grid, the regions are from HadCM2 

subroutine SelectRS

do
  call GridSelect (RSGridChosen,RSGridTitle,RSGridLongN,RSGridLatN,RSGridDataN,RSGridFile)
  if (RSGridLongN.NE.ExeN.OR.RSGridLatN.NE.WyeN) &
  	print "(a,4i5)", "   > Selected grid does not match. Re-try.", RSGridLongN,ExeN,RSGridLatN,WyeN
  if (RSGridLongN.EQ.ExeN.AND.RSGridLatN.EQ.WyeN) exit
end do

call RegSelect (RSGridChosen, RSGridLongN, RSGridLatN, RSGridDataN, RSMapIDLReg, RSRegSizes, RSRegNames, &
			RSRegTitle, RSRegN, RegFileOut=CommInfo(XComm,1))    

deallocate (RSMapIDLReg, RSRegSizes, RSRegNames, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SelectRS: Deallocation failure #####"

end subroutine SelectRS

!*******************************************************************************
! select truncation limits --> 4,5

subroutine SelectTL

print*, "  > Specify the lower, upper limits (-999=no limit):"
do
  read (*,*,iostat=ReadStatus), CommOp(XComm,4),CommOp(XComm,5)
  if (CommOp(XComm,4).NE.MissVal.AND.CommOp(XComm,5).NE.MissVal.AND.CommOp(XComm,4).GT.CommOp(XComm,5)) then
	  	ReadStatus = 1
	  	print*, "  > Unacceptable limits. Try again."
  end if
  if (ReadStatus.LE.0) exit
end do
    
end subroutine SelectTL

!*******************************************************************************
! decide whether to weight normally (=1,each box=1) or by area (=2)
! must follow SelectRS, because needs CommInfo(XComm,1) to be filled

subroutine SelectWE

do
	read (*,*,iostat=ReadStatus), QWeight
	if (ReadStatus.LE.0.AND.QWeight.GE.1.AND.QWeight.LE.2) exit
end do
CommOp(XComm,4) = QWeight

if (QWeight.EQ.2) then			! ensures that specified weights file actually exists
  GivenFile  = GetGloRef (ExeN,WyeN)
  OrigSubBeg = index(GivenFile,".ref")
  GivenFile  = GivenFile(1:(OrigSubBeg-1)) // ".weights.glo"
  print*, "  > You must ensure that a weights .glo exists named: "	! the .glo must be rpb, poles incl
  print "(2a)", "   > ", trim(GivenFile)
end if

end subroutine SelectWE

!*******************************************************************************
! specify automatics

subroutine AutoSpec

print*, "  > Enter the substring in the filepaths to vary: "
do
	read (*,*,iostat=ReadStatus), OrigSub
	if (ReadStatus.GT.0) then	
			print*, "  > Bad format. Try again."
	else if (OrigSub.EQ."") then
			print*, "  > Blank not permitted. Try again."
	end if
	if (ReadStatus.LE.0.AND.OrigSub.NE."") exit
end do
	
SubLen = len(trim(OrigSub))
  
print*, "  > Enter the substring in each execution, starting with no.2:"
do XExec = 2, ExecN
  do
    read (*,*,iostat=ReadStatus), SpecSub
    if (ReadStatus.GT.0) then	
			print*, "  > Bad format. Try again."
    else if (SpecSub.EQ."") then
			print*, "  > Blank not permitted. Try again."
    end if
    if (ReadStatus.LE.0.AND.SpecSub.NE."") exit
  end do
        
  do XComm = 1, CommN
    if (CommFile(XComm,1).NE."") then						! provide file
      GivenFile  = CommFile(XComm,1)
      FileSubBeg = index(GivenFile,OrigSub(1:SubLen))				! start of orig sub
										! replace orig sub if poss
      if (FileSubBeg.GT.0) CommFile(XComm,XExec) = GivenFile(1:(FileSubBeg-1)) // trim(SpecSub) // &
      					trim(GivenFile((FileSubBeg+SubLen):80))
            
      if (CommInfo(XComm,1).NE."") then						! provide info line
              CommInfo(XComm,XExec) = CommInfo(XComm,1)
      end if
    end if
  end do
end do
  
end subroutine AutoSpec

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

subroutine SaveSpecFile

print*, "  > Saving to .ops ..."

open  (1,file=SpecFile,status="new",access="sequential",form="formatted",action="write")

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

do XComm = 1, CommN
  write (1,"(5i9,f9.2)"), (CommOp(XComm,XElement), XElement=1,5), CommKay(XComm)

  do XExec = 1, ExecN
    write (1,"(a80)"), CommFile(XComm,XExec)
    write (1,"(a80)"), CommInfo(XComm,XExec)
  end do
end do

do XExe = 1, ExeN
  do XWye = 1, WyeN
    write (1,"(i10)"), RefGrid(XExe,XWye)
  end do
end do

close (1)

end subroutine SaveSpecFile

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

subroutine Finalise

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

print*

end subroutine Finalise

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

end program SetupSeqGrim