! tynsc20seqgrim_pik.f90
! f90 -o tynsc20seqgrim_pik time_pik.f90 filenames_pik.f90 grimfiles_pik.f90 annfiles_pik.f90 basicfun_pik.f90 tynsc20seqgrim_pik.f90
! written by Tim Mitchell on 30.03.01 and subsequently modified
! it is designed to be run using Comms defined using a specification file from *unpack*.f90
! This file has been modified by Soenke Zaehle (PIK) on 30.07.02 for use with XLF compiler
! (PIK). changes are: 
! - removed "," after read and write statements
! - replaced SCRATCH statements with REPLACE statements

program TYNsc20SeqGrim

use Time
use FileNames
use GrimFiles
use AnnFiles
use BasicFun

implicit none

real, pointer, dimension (:,:,:)		:: OperA, OperB, OperC, OperD
real, pointer, dimension (:,:,:)		:: Grim1, Grim2, Grim3, Grim4, Grim5, Grim6
real, pointer, dimension (:,:,:)		:: GrimXtra
real, pointer, dimension (:,:)			:: GridWeights, WeightsXtra, LinData,MonthData,SeasonData
real, pointer, dimension (:)			:: CommKay, VecA, VecB, VecC, VecD, GloSlice,AnnualData
real, pointer, dimension (:)			:: RegWeights,GloWeights
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, MonthLengths
integer, pointer, dimension (:)			:: YearAD, FileYearAD

character (len=80), pointer, dimension (:,:)	:: CommFile, CommInfo
character (len=80), pointer, dimension (:)	:: FileLineNames,RegFilePaths,OpsFile
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    = ""

real :: OpTot, OpEn, OpMiss, OpMean, OpThresh, OpValid, OpBelow, OpAbove
real :: FileFactor, LSRAye,LSRBee,LSRPea

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
integer :: PerYear0, PerYear1, YearAD0, YearAD1, FileYear0, FileYear1, MasterYear0, MasterYear1
integer :: ThisMonth, ThisYear, NoZip, TotA,TotB,TotC,TotD, OpsFileN,Operate
integer :: QMeanSum

character (len=80) :: FileInfo,FilePath,FileVariable,GridFilePath,GloFile,GloInfo,GenericFile,SaveVariName
character (len=80) :: Stem,Tip,ExtraStr,SpecFile,Filter,GridFile
character (len=20) :: RegName
character (len=40) :: RegTitle
character (len= 4) :: FileSuffix, SaveSuffix

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

open  (99,file="log.dat",status="replace",action="write")

call FindLatestOpsFile
call Initialise
call Main
call Finalise

close (99)

contains

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

subroutine FindLatestOpsFile

Filter="*.ops"
call GetBatch (Filter,OpsFile)
OpsFileN = size(OpsFile)
SpecFile = OpsFile(OpsFileN)

end subroutine FindLatestOpsFile

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

subroutine Initialise

nullify (OperA,OperB,OperC,OperD)
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 "(2a)", "   > Operating on: ", trim(SpecFile)
print* 

call LoadSpec

allocate (Grim1 (YearN,MonthN,BoxN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Initialise: Allocation failure: Grim1 #####"
  
allocate (Grim2 (YearN,MonthN,BoxN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: Initialise: Allocation failure: Grim2 #####"

end subroutine Initialise

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

subroutine LoadSpec

open  (1,file=SpecFile,status="old",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,*) "SeqGrim : data found within comms file."
do XComm = 1, CommN
  read (1,"(5i9,f9.2)")  (CommOp(XComm,XElement), XElement=1,5), CommKay(XComm)
  write (99,"(6i6,f12.4)") XComm, (CommOp(XComm,XElement), XElement=1,5), CommKay(XComm)
  
  do XExec = 1, ExecN
    read (1,"(a)") CommFile(XComm,XExec)
    read (1,"(a)") CommInfo(XComm,XExec)
    write (99,"(a)") trim(adjustl(CommFile(XComm,XExec)))
    write (99,"(a)") trim(adjustl(CommInfo(XComm,XExec)))
  end do
end do
write (99,*) "SeqGrim : data found within comms file. END."

close (1)

allocate (YearAD (YearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadSpec: Allocation failure: G+YAD #####"
do XYear = 1, YearN
  YearAD(XYear) = YearAD0 + XYear - 1
end do

GridFile = 'obs.clim6190.ateam.cld'
call LoadGrim (GrimXtra,Grid,FileYearAD,FileBounds,FileInfo,GridFile,'.cld',FileSuffix)

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

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. 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.10) then
        print "(a,i1,a,i4,a1,i4,a,i1)", "   > Climatology: ", CommOp(XComm,2), " = ", CommOp(XComm,4), "-", &
        			CommOp(XComm,5), " mean from ", CommOp(XComm,3)
    	call OpClim
    else if (CommOp(XComm,1).EQ.17) then
    	call OpMaskRel
        print "(a,i1,a,i2,a,i1,a,i1,a,i1,a,i7)", "   > 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.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.43) then
        print "(a,i1)", "   > Saving from: ", CommOp(XComm,2)
	call DumpGrim
    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), &
			": success=", TotA, ", failure=", TotB
    end if
  end do
end do

end subroutine Main

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

subroutine ResetAll

if (associated(Grim1)) Grim1=MissVal
if (associated(Grim2)) Grim2=MissVal

end subroutine ResetAll

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

subroutine FindVariable

XComm = 0 ; FileSuffix = "" ; FileVariable = "" ; FileFactor = MissVal
do
  XComm = XComm + 1

  if (CommFile(XComm,XExec).NE."") then
  	SuffixStart = index(CommFile(XComm,XExec),'.',.TRUE.)
  	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 (FileSuffix.EQ.".pre".OR.FileSuffix.EQ.".frs".OR.FileSuffix.EQ.".wet") then
  QMeanSum = 2
else
  QMeanSum = 0
end if

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

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
! this has been altered to ignore the year structure and make it 2001-end

subroutine GrabGrim

call PointOperA
print "(2a)", "   > Seeking: ", trim(CommFile(XComm,XExec))
call LoadGrim (OperA,FileGrid,FileYearAD,FileBounds,FileInfo,CommFile(XComm,XExec),"    ",FileSuffix,&
		MasterYearAD=YearAD,OverRide=2001)

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)) CheckGrid = CheckGrid + 1
  end do
end do
if (CheckGrid.NE.0) print*, "  > ##### ERROR: Grid mismatches: ", CheckGrid

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
print "(2a)", "   > Seeking: ", trim(CommFile(XComm,XExec))
call LoadGrim (GrimXtra,FileGrid,FileYearAD,FileBounds,FileInfo,CommFile(XComm,XExec),"    ",FileSuffix)

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)) CheckGrid = CheckGrid + 1
  end do
end do
if (CheckGrid.NE.0) print*, "  > ##### ERROR: Grid mismatches: ", CheckGrid

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 = 3

subroutine GrabAnn

call PointOperA
print "(2a)", "   > Seeking: ", trim(CommFile(XComm,XExec))
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
 
if (PerYear1.GT.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 (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 = 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 = 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 = 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 = 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

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

subroutine PointOperA

if (CommOp(XComm,2).EQ.1) OperA => Grim1
if (CommOp(XComm,2).EQ.2) OperA => Grim2

end subroutine PointOperA

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

subroutine PointOperB

if (CommOp(XComm,3).EQ.1) OperB => Grim1
if (CommOp(XComm,3).EQ.2) OperB => Grim2

end subroutine PointOperB

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

subroutine PointOperC

if (CommOp(XComm,5).EQ.1) OperC => Grim1
if (CommOp(XComm,5).EQ.2) OperC => Grim2

end subroutine PointOperC

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

subroutine PointOperD

if (CommOp(XComm,4).EQ.1) OperD => Grim1
if (CommOp(XComm,4).EQ.2) OperD => Grim2

end subroutine PointOperD

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

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 #####"

print*

end subroutine Finalise

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

end program TYNsc20SeqGrim