! globulk.f90
! performs basic operations on .glo files in bulk
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
! 	-o ./../goglo/globulk filenames.f90 sortmod.f90 initialmod.f90 glofiles.f90 
!    	gloops.f90 basicfun.f90 ./../goglo/globulk.f90 2> /tyn1/tim/scratch/stderr.txt
! last modified 10.01.02

program GloBulk

use FileNames
use SortMod
use InitialMod
use GloFiles
use GloOps
use BasicFun

implicit none

!*******************************************************************************
! definitions

real, pointer, dimension (:)		:: RegVector,MidVector,OutVector,RegWeights
real, pointer, dimension (:,:)		:: GridWeights,FileRegData,Deg025,Deg050,DegNew

integer, pointer, dimension (:,:) 	:: MapIDLReg,MintMapIDLReg,OrigMapIDLReg  ! shape of IDL mapping arrays
integer, pointer, dimension (:,:) 	:: Biome
integer, pointer, dimension (:) 	:: RegSizes, MintRegSizes,OrigRegSizes,Order
integer, dimension (16) 		:: BiomeReg

character (len=20), pointer, dimension (:) :: RegNames,MintRegNames,OrigRegNames  ! names of individual regions
character (len=80), pointer, dimension (:) :: InFile, MidFile, OutFile

real, parameter :: MissVal = -999.0

real :: OpTot, OpMiss, OpEn, OpConstant, HugeVal, TinyVal, MagMin,MagMax, OpVal
real :: Lower,Upper,OpMin,OpMax

integer :: ReadStatus, AllocStat
integer :: FileN,RegN,MidFileN
integer :: XFile,XReg,XMidFile,XLong,XLat,XMiniLong,XMiniLongAny,XMiniLat,XPre,XTmp
integer :: XLat050,XLon050,XLat025,XLon025,XOut050,XLatRJ,XLonRJ,XLatAve,XLonAve,XExe,XWye
integer :: QOperation, QFinish, QOpPerform, QWeight, QStat, QOutMag, QReplicate
integer :: InSubLen,InSubBeg,InFileLen,LastSlash,LastMiss,MiniRad,NUnique,NMiss
integer :: GridChosen,GridLongN,GridLatN,GridDataN, BiomeOther,BiomeOut			

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

character (len=80) :: GivenFile, NameOnly, GridFile, GloSaveFile, BiomeFile
character (len=80) :: InSub, OutSub, SaveTitle
character (len=40) :: RegTitle			
character (len=10) :: GridTitle			

!*******************************************************************************
! main program

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

print*
print*, "  > ***** GloBulk: performs basic ops on .glo files in bulk *****"
print*

call SelectOp

print*

close (99)

contains

!*******************************************************************************
! select operation to perform on each .glo

subroutine SelectOp

print*, "  > Select the operation to perform (0=list): "
QOperation = -1
QFinish    = -1

do
  if      (QOperation.EQ.0) then
    print*, "  >  1. transform into different region set (weighting=option)"
    print*, "  >  2. place supra-region mean in dump file" 
    print*, "  >  3. fill gaps and save to new .glo" 
    print*, "  >  4. place supra-file median in .glo" 
    print*, "  >  5. place supra-file mean   in .glo" 
    print*, "  >  6. place supra-file count  in .glo" 
    print*, "  >  7. place supra-file classd in .glo" 
    print*, "  >  8. operate: .glo . constant"
    print*, "  >  9. operate: .glo . .glo"
    print*, "  > 11. convert RichJ to 0.5deg"
    print*, "  > 12. convert climate to biomes"
    print*, "  > 13. report grid-box value"
    QOperation = -1
  else if (QOperation.EQ.1) then
    call NewRegSet ; QFinish = 1
  else if (QOperation.EQ.2) then
    call SupraRegMean ; QFinish = 1
  else if (QOperation.EQ.3) then
    call FillGapsSave ; QFinish = 1
  else if (QOperation.EQ.4) then
    QStat = 1 ; call SupraFileMed ; QFinish = 1
  else if (QOperation.EQ.5) then
    QStat = 2 ; call SupraFileMed ; QFinish = 1
  else if (QOperation.EQ.6) then
    QStat = 3 ; call SupraFileMed ; QFinish = 1
  else if (QOperation.EQ.7) then
    QStat = 4 ; call SupraFileMed ; QFinish = 1
  else if (QOperation.EQ.8) then
    call OpGloConstant ; QFinish = 1
  else if (QOperation.EQ.9) then
    call OpMultiGlo ; QFinish = 1
  else if (QOperation.EQ.11) then
    call OpRichJHalfDeg ; QFinish = 1
  else if (QOperation.EQ.12) then
    call ConvertBiome ; QFinish = 1
  else if (QOperation.EQ.13) then
    call ReportGridBox ; QFinish = 1
  else
    do
      read (*,*,iostat=ReadStatus), QOperation
      if (ReadStatus.LE.0) exit
    end do
  end if
  
  if (QFinish.EQ.1) exit
end do

end subroutine SelectOp

!*******************************************************************************
!  1. transform into different region set

subroutine NewRegSet

print*, "  > Identify the OLD region set (constants/goglo/ref/)."
call LoadRefReg (Blank,OrigMapIDLReg,OrigRegSizes,OrigRegNames)

print*, "  > Identify the NEW region set (constants/goglo/ref/)."
call LoadRefReg (Blank,MintMapIDLReg,MintRegSizes,MintRegNames)

print*, "  > Weight the aggregations by box (=1) or area (=2) ?"
do
      read (*,*,iostat=ReadStatus), QWeight
      if (ReadStatus.LE.0.AND.QWeight.GE.1.AND.QWeight.LE.2) exit
end do

GridLongN = size(OrigMapIDLReg,1) ; GridLatN = size(OrigMapIDLReg,2)

if (QWeight.EQ.1) then
  allocate (GridWeights(GridLongN,GridLatN),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: NewRegSet: Allocation failure: Weights #####"
  GridWeights = 1.0
else
  GivenFile  = GetGloRef (GridLongN,GridLatN)
  InSubBeg   = index(GivenFile,".ref")
  GivenFile  = GivenFile(1:(InSubBeg-1)) // ".weights.glo"
  			! what we want is a weights .glo of rpb form, incl polar boxes 
  print "(2a)", "   > We look for: ", trim(adjustl(GivenFile))
  call LoadGloGrid (GivenFile, GridWeights, Silent=1)
  if (size(GridWeights,1).NE.GridLongN.OR.size(GridWeights,2).NE.GridLatN) then
    GridWeights = MissVal
    print*, "  > ##### ERROR: NewRegSet: Loaded .glo has wrong dims #####"
  end if
end if

call CheckOutMag

print*, "  > Select the OLD files. There must be a common substring."
call GetBatch (Blank,InFile)
FileN  = size (InFile, 1)

call GetOutFile

do XFile = 1, FileN
  GivenFile = InFile(XFile)
  LastSlash = index(GivenFile,"/",.TRUE.)
  NameOnly  = adjustl(trim(GivenFile((LastSlash+1):80)))
  print "(a,i4,2a)", "   > Batch item: ", XFile, ": ", trim(NameOnly)
  if (QOutMag.EQ.0) then
    call GloToGlo (InFile(XFile),OutFile(XFile),OrigMapIDLReg,MintMapIDLReg,CallWeights=GridWeights)
  else
    call GloToGlo (InFile(XFile),OutFile(XFile),OrigMapIDLReg,MintMapIDLReg,CallWeights=GridWeights, &
    			MagMin=MagMin,MagMax=MagMax)
  end if
end do

deallocate (OrigMapIDLReg,OrigRegSizes,OrigRegNames, &
	    MintMapIDLReg,MintRegSizes,MintRegNames, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: NewRegSet: Deallocation failure: most #####"

deallocate (InFile,OutFile,GridWeights,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: NewRegSet: Deallocation failure: file #####"

end subroutine NewRegSet

!*******************************************************************************
!  2. put supra-regional mean in dump file

subroutine SupraRegMean

print*, "  > Specify the common grid. "
call GridSelect (GridChosen,GridTitle,GridLongN,GridLatN,GridDataN,GridFile)
call RegSelect (GridChosen,GridLongN,GridLatN,GridDataN,MapIDLReg, RegSizes, RegNames, RegTitle, RegN)

print*, "  > Select the OLD files. There must be a common substring."
call GetBatch (Blank,InFile)
FileN  = size (InFile, 1)

open (98,file="./../../../scratch/dump-globulk.dat",status="replace",action="write")

do XFile = 1, FileN
  GivenFile = InFile(XFile)
  LastSlash = index(GivenFile,"/",.TRUE.)
  NameOnly  = adjustl(trim(GivenFile((LastSlash+1):80)))
  print "(a,i4,2a)", "   > Batch item: ", XFile, ": ", trim(NameOnly)
  
  call LoadGloVec (InFile(XFile), MapIDLReg, RegVector, Silent=1)
  
  OpTot=0.0 ; OpMiss=0.0 ; OpEn=0.0
  do XLong = 1, GridLongN
    do XLat = 1, GridLatN
      if (MapIDLReg(XLong,XLat).NE.MissVal) then
        OpEn=OpEn+1
        if (RegVector(MapIDLReg(XLong,XLat)).NE.MissVal) then
          OpTot=OpTot+RegVector(MapIDLReg(XLong,XLat)) 
        else
          OpMiss=OpMiss+1
        end if
      end if
    end do
  end do
  
  if (OpEn.GT.0.AND.(OpEn-OpMiss).GT.0) then
    write (98,"(2e12.5,a80)"), (OpTot/(OpEn-OpMiss)), (OpMiss/OpEn), trim(NameOnly)
  else
    write (98,"(2e12.5,a80)"), MissVal, MissVal, trim(NameOnly)
  end if  
  
  deallocate (RegVector, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: SupraRegMean: Deallocation failure: RV #####"
end do

close (98)

deallocate (MapIDLReg, RegSizes, RegNames, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SupraRegMean: Deallocation failure: all #####"

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

end subroutine SupraRegMean

!*******************************************************************************
!  3. fill gaps in file and resave

subroutine FillGapsSave

print*, "  > Specify the common grid. "
call GridSelect (GridChosen,GridTitle,GridLongN,GridLatN,GridDataN,GridFile)
call RegSelect (GridChosen,GridLongN,GridLatN,GridDataN,MapIDLReg, RegSizes, RegNames, RegTitle, RegN)

print*, "  > Select the OLD files. There must be a common substring."
call GetBatch (Blank,InFile)
FileN  = size (InFile, 1)

call GetOutFile

do XFile = 1, FileN
  GivenFile = InFile(XFile)
  LastSlash = index(GivenFile,"/",.TRUE.)
  NameOnly  = adjustl(trim(GivenFile((LastSlash+1):80)))
  print "(a,i4,2a)", "   > Batch item: ", XFile, ": ", trim(NameOnly)
  
  call LoadGloVec (InFile(XFile), MapIDLReg, RegVector, Silent=1)
  
  LastMiss = GridLongN * GridLatN ; MiniRad = 1
  do 
   OpMiss = 0
   
   do XLong = 1, GridLongN
    do XLat = 1, GridLatN
      if (MapIDLReg(XLong,XLat).NE.MissVal) then
        if (RegVector(MapIDLReg(XLong,XLat)).EQ.MissVal) then
          OpTot=0.0 ; OpEn=0.0
          
          do XMiniLongAny = (XLong-MiniRad), (XLong+MiniRad)
            XMiniLong = XMiniLongAny
            if (XMiniLong.LT.1) XMiniLong = GridLongN + XMiniLong
            if (XMiniLong.GT.GridLongN) XMiniLong = XMiniLong - GridLongN 
            
            do XMiniLat = (XLat-MiniRad), (XLat+MiniRad)
             if (XMiniLat.GE.1.AND.XMiniLat.LE.GridLatN) then
              if (MapIDLReg(XMiniLong,XMiniLat).NE.MissVal) then
                if (RegVector(MapIDLReg(XMiniLong,XMiniLat)).NE.MissVal) then
                  OpTot = OpTot + RegVector(MapIDLReg(XMiniLong,XMiniLat))
                  OpEn  = OpEn  + 1
                end if
              end if
             end if
            end do
          end do
          if (OpEn.GE.1) then
          	RegVector(MapIDLReg(XLong,XLat)) = OpTot / OpEn
          else
          	OpMiss = OpMiss + 1
          end if
        end if
      end if
    end do
   end do
   
   if (OpMiss.EQ.LastMiss) MiniRad = MiniRad + 1
   LastMiss = OpMiss
   
   if (OpMiss.EQ.0) exit
  end do
  
  SaveTitle = trim(NameOnly) // " : " // "with gaps filled"
  call SaveGlo (GridLongN,GridLatN,RegN,GridFile,OutFile(XFile),SaveTitle,RegVector,MapIDLReg)
  
  deallocate (RegVector, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: FillGapsSave: Deallocation failure: RV #####"
end do

deallocate (MapIDLReg, RegSizes, RegNames, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: FillGapsSave: Deallocation failure: all #####"

deallocate (InFile,OutFile,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: FillGapsSave: Deallocation failure: file #####"

end subroutine FillGapsSave

!*******************************************************************************
!  4,5,6,7 place supra-file median,mean,count,classd in .glo
! count gives the number of unique values for a region, from the set of files
! classd gives a .glo containing the number of .glo with values in a certain range
!   e.g. 4,4,8,1 would give 3

subroutine SupraFileMed

print*, "  > Specify the common grid. "
call GridSelect (GridChosen,GridTitle,GridLongN,GridLatN,GridDataN,GridFile)
call RegSelect (GridChosen,GridLongN,GridLatN,GridDataN,MapIDLReg, RegSizes, RegNames, RegTitle, RegN)

print*, "  > Select the OLD files. There must be a common substring."
call GetBatch (Blank,InFile)
FileN  = size (InFile, 1)

if (QStat.EQ.4) then
  print*, "  > Enter the lower,upper bounds (-999=none):"
  do
      read (*,*,iostat=ReadStatus), Lower,Upper
      if (ReadStatus.LE.0.AND.(Lower.EQ.MissVal.OR.Upper.EQ.MissVal.OR.Lower.LE.Upper)) exit
  end do
end if

allocate (FileRegData(FileN,RegN),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SupraFileMed: Allocation failure: FRD #####"

do XFile = 1, FileN
  GivenFile = InFile(XFile)
  LastSlash = index(GivenFile,"/",.TRUE.)
  NameOnly  = adjustl(trim(GivenFile((LastSlash+1):80)))
  print "(a,i4,2a)", "   > Batch item: ", XFile, ": ", trim(NameOnly)
  
  call LoadGloVec (InFile(XFile), MapIDLReg, RegVector, Silent=1)
    
  do XReg = 1, RegN
    FileRegData(XFile,XReg) = RegVector(XReg)
  end do
  
  deallocate (RegVector,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: SupraFileMed: Deallocation failure: RV #####"
end do

allocate (MidVector(FileN),&
	  OutVector(RegN),Order(FileN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SupraFileMed: Allocation failure: MV #####"

if (QStat.EQ.1.OR.QStat.EQ.2) then
 do XReg = 1, RegN
  do XFile = 1, FileN
    MidVector(XFile) = FileRegData(XFile,XReg)
  end do
  
  if (QStat.EQ.1) then
    call QuickSort (Reals=MidVector,OrderValid=Order,NMiss=NMiss)
    OutVector(XReg) = MidVector (Order (nint((FileN-NMiss)/2.0)) )
  else if (QStat.EQ.2) then
    OutVector(XReg) = OpCalcMean  (MidVector,FileN,10.0)
  end if
 end do
else if (QStat.EQ.3) then
  do XReg = 1, RegN
    MidVector=MissVal ; NUnique=0
    do XFile = 1, FileN
      XMidFile=0 
      do
        XMidFile=XMidFile+1
	if (FileRegData(XFile,XReg).NE.MissVal.AND.MidVector(XMidFile).EQ.MissVal) then
	  MidVector(XMidFile)=FileRegData(XFile,XReg)
	  NUnique=NUnique+1
	end if
	if (MidVector(XMidFile).NE.MissVal.AND. &
		MidVector(XMidFile).EQ.FileRegData(XFile,XReg)) exit
	if (XMidFile.EQ.FileN) exit
      end do
    end do
    OutVector(XReg) = NUnique
  end do 
else if (QStat.EQ.4) then
 OutVector=0
 do XReg = 1, RegN
  do XFile = 1, FileN
    if (FileRegData(XFile,XReg).NE.MissVal) then
      if (Lower.EQ.MissVal.OR.FileRegData(XFile,XReg).GE.Lower) then
        if (Upper.EQ.MissVal.OR.FileRegData(XFile,XReg).LE.Upper) then
	  OutVector(XReg)=OutVector(XReg)+1
	end if
      end if
    end if
  end do
 end do
end if

OpMax=maxval(outvector) ; OpMin=minval(outvector)

print*, "  > Enter the .glo file to save (range:",OpMin,OpMax
do
      read (*,*,iostat=ReadStatus), GivenFile
      if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
GloSaveFile = SavePath (GivenFile,".glo")

print*, "  > Enter the information line:"
do
      read (*,*,iostat=ReadStatus), SaveTitle
      if (ReadStatus.LE.0.AND.SaveTitle.NE."") exit
end do

call SaveGlo (GridLongN,GridLatN,RegN,GridFile,GloSaveFile,SaveTitle,OutVector,MapIDLReg)  

deallocate (FileRegData,MidVector,OutVector,MapIDLReg,RegSizes,RegNames,Order,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SupraFileMed: Deallocation failure: RV #####"

end subroutine SupraFileMed

!*******************************************************************************
!  8. .glo . constant

subroutine OpGloConstant

print*, "  > Select the op: A/k=1,A*k=2,A-k=3,A+k=4,sqroot(A)=5,A**k=6,abs(A)=7"
do
      read (*,*,iostat=ReadStatus), QOpPerform
      if (ReadStatus.LE.0.AND.QOpPerform.GE.1.AND.QOpPerform.LE.7) exit
end do

print*, "  > Select the constant:"
do
      read (*,*,iostat=ReadStatus), OpConstant
      if (ReadStatus.LE.0) exit
end do

call CheckOutMag

print*, "  > Specify the common grid. "
call GridSelect (GridChosen,GridTitle,GridLongN,GridLatN,GridDataN,GridFile)
call RegSelect (GridChosen,GridLongN,GridLatN,GridDataN,MapIDLReg, RegSizes, RegNames, RegTitle, RegN)

print*, "  > Select the OLD files. There must be a common substring."
call GetBatch (Blank,InFile)
FileN  = size (InFile, 1)

call GetOutFile

do XFile = 1, FileN
  GivenFile = InFile(XFile)
  LastSlash = index(GivenFile,"/",.TRUE.)
  NameOnly  = adjustl(trim(GivenFile((LastSlash+1):80)))
  print "(a,i4,2a)", "   > Batch item: ", XFile, ": ", trim(NameOnly)
  
  call LoadGloVec (InFile (XFile), MapIDLReg, RegVector, Silent=1)
  allocate (OutVector(RegN),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: OpGloConstant: Allocation failure: OV #####"
  
  do XReg = 1, RegN
    OutVector(XReg) = OpAwithB (RegVector(XReg),OpConstant,QOpPerform)
  end do
  
  if (QOutMag.EQ.1) then
   do XReg = 1, RegN
    if (OutVector(XReg).NE.MissVal) then
      if (abs(OutVector(XReg)).LT.MagMin) OutVector(XReg) = sign(MagMin,OutVector(XReg))
      if (abs(OutVector(XReg)).GT.MagMax) OutVector(XReg) = sign(MagMax,OutVector(XReg))
    end if
   end do
  end if
  
  SaveTitle = trim(NameOnly) // " operated on"
  call SaveGlo (GridLongN,GridLatN,RegN,GridFile,OutFile(XFile),SaveTitle,OutVector,MapIDLReg)
  
  deallocate (RegVector,OutVector, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: OpGloConstant: Deallocation failure: RV #####"
end do  

deallocate (MapIDLReg, RegSizes, RegNames, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: OpGloConstant: Deallocation failure: all #####"

deallocate (InFile,OutFile,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: OpGloConstant: Deallocation failure: file #####"

end subroutine OpGloConstant

!*******************************************************************************
!  9. operate one .glo on another

subroutine OpMultiGlo

QReplicate = 0

print*, "  > Select the op: A/B=1,A*B=2,A-B=3,A+B=4,sqroot(A)=5,A**B=6,abs(A)=7"
do
      read (*,*,iostat=ReadStatus), QOpPerform
      if (ReadStatus.LE.0.AND.QOpPerform.GE.1.AND.QOpPerform.LE.7) exit
end do

call CheckOutMag

print*, "  > Specify the common grid. "
call GridSelect (GridChosen,GridTitle,GridLongN,GridLatN,GridDataN,GridFile)
call RegSelect (GridChosen,GridLongN,GridLatN,GridDataN,MapIDLReg, RegSizes, RegNames, RegTitle, RegN)
call GetParallelFiles
call GetOutFile

do XFile = 1, FileN
  call GetMidFileIndex
  
  call LoadGloVec (InFile (XFile), MapIDLReg, RegVector, Silent=1)	! load 'in' file
  call LoadGloVec (MidFile(XMidFile), MapIDLReg, MidVector, Silent=1)  	! load 'mid' file
  
  allocate (OutVector(RegN),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: OpMultiGlo: Allocation failure: OV #####"
  
  do XReg = 1, RegN
    OutVector(XReg) = OpAwithB (RegVector(XReg),MidVector(XReg),QOpPerform)
  end do
  
  if (QOutMag.EQ.1) then
   do XReg = 1, RegN
    if (OutVector(XReg).NE.MissVal) then
      if (abs(OutVector(XReg)).LT.MagMin) OutVector(XReg) = sign(MagMin,OutVector(XReg))
      if (abs(OutVector(XReg)).GT.MagMax) OutVector(XReg) = sign(MagMax,OutVector(XReg))
    end if
   end do
  end if
  
  SaveTitle = 'result of bulk operation'
  call SaveGlo (GridLongN,GridLatN,RegN,GridFile,OutFile(XFile),SaveTitle,OutVector,MapIDLReg)
  
  deallocate (RegVector,MidVector,OutVector, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: OpMultiGlo: Deallocation failure: RV #####"
end do  

deallocate (MapIDLReg, RegSizes, RegNames,InFile,OutFile,MidFile,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: OpMultiGlo: Deallocation failure: file #####"

end subroutine OpMultiGlo

!*******************************************************************************
! 11. convert RichJ .glo to 0.5deg .glo

subroutine OpRichJHalfDeg

print*, "  > Select OLD files on RichJ grid. There must be a common substring."
call GetBatch (Blank,InFile)
FileN  = size (InFile, 1)

call GetOutFile

GridFile = "./../../../constants/goglo/ref/half-degree.ref"

allocate (Deg025(720,1440), &
	  DegNew(720,1440), &
	  Deg050(360, 720), &
	  OutVector(259200),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: OpRichJHalfDeg: Allocation failure #####"

do XFile = 1, FileN
  GivenFile = OutFile(XFile)
  LastSlash = index(GivenFile,"/",.TRUE.)
  NameOnly  = adjustl(trim(GivenFile((LastSlash+1):80)))
  print "(a,i4,2a)", "   > Batch item: ", XFile, ": ", trim(NameOnly)
  
  call LoadGloGrid (InFile (XFile), FileRegData, Silent=1)
  
  allocate (MapIDLReg(720,360), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: OpRichJHalfDeg: Reallocation failure: MapIDLReg #####"

  Deg025=MissVal ; Deg050=MissVal ; DegNew=MissVal
  
  XLatRJ=1
  do XLat025=1,720
    XLonRJ=1
    if (mod((real(XLat025)-6.0),10.0).EQ.0.0) XLatRJ=XLatRJ+1
    do XLon025=1,1440
      if (mod((real(XLon025)-0.0),10.0).EQ.0.0) XLonRJ=XLonRJ+1      
      Deg025(XLat025,XLon025) = FileRegData(XLonRJ,XLatRJ)
    end do
  end do
  
  do XLat025=1,720
    do XLon025=1,1440
      OpTot=0.0 ; OpEn=0.0
      do XLatAve=max(1,XLat025-4),min(720,XLat025+4)
        do XLonAve=XLon025-4,XLon025+4
          if (XLonAve.LT.1) then
            if (Deg025(XLatAve,XLonAve+1440).NE.MissVal) then
              OpTot=OpTot+Deg025(XLatAve,XLonAve+1440) ; OpEn=OpEn+1
            end if
          else if (XLonAve.GT.1440) then
            if (Deg025(XLatAve,XLonAve-1440).NE.MissVal) then
              OpTot=OpTot+Deg025(XLatAve,XLonAve-1440) ; OpEn=OpEn+1
            end if
          else
            if (Deg025(XLatAve,XLonAve).NE.MissVal) then
              OpTot=OpTot+Deg025(XLatAve,XLonAve) ; OpEn=OpEn+1
            end if
          end if
        end do        
      end do
      if (OpEn.GE.09) DegNew(XLat025,XLon025) = OpTot/OpEn
    end do
  end do
  
  do XLat050=1,360
    do XLon050=1,720
      OpTot=0.0 ; OpEn=0.0
      do XLat025=((XLat050-1)*2)+1,((XLat050-1)*2)+2
        do XLon025=((XLon050-1)*2)+1,((XLon050-1)*2)+2
          if (DegNew(XLat025,XLon025).NE.MissVal) then
            OpTot=OpTot+DegNew(XLat025,XLon025) ; OpEn=OpEn+1
          end if
        end do
      end do
      if (OpEn.GT.0) Deg050(XLat050,XLon050) = OpTot/OpEn
    end do
  end do
  
  XOut050=0
  do XLon050=1,720
    do XLat050=1,360
      XOut050=XOut050+1
      MapIDLReg(XLon050,XLat050)=XOut050
      OutVector(XOut050)=Deg050(XLat050,XLon050)
    end do
  end do

  SaveTitle = '144*73->720*360'
  call SaveGlo (720,360,259200,GridFile,OutFile(XFile),SaveTitle,OutVector,MapIDLReg)
  
  deallocate (FileRegData,MapIDLReg, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: OpRichJHalfDeg: Deallocation failure: loop #####"
end do

deallocate (Deg025,Deg050,DegNew,OutVector,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: OpRichJHalfDeg: Deallocation failure #####"

end subroutine OpRichJHalfDeg

!*******************************************************************************
! 12. convert climate to biomes

subroutine ConvertBiome

QReplicate = 0 

allocate (Biome(45,45),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertBiome: Allocation failure: OV #####"
Biome=MissVal
BiomeFile="./../../../constants/biome/biomedefine.txt"

open  (2,file=BiomeFile,status="old",access="sequential",form="formatted",action="read")
do XPre=1,17
  read (2, fmt="(a)"), GivenFile	! just a convenient trash form
end do
do XTmp=45,1,-1
  read (2, fmt="(45i3)"), (Biome(XTmp,XPre),XPre=1,45) 
end do
close (2)

print*, "  > Specify the common grid. "
call GridSelect (GridChosen,GridTitle,GridLongN,GridLatN,GridDataN,GridFile)
call RegSelect (GridChosen,GridLongN,GridLatN,GridDataN,MapIDLReg, RegSizes, RegNames, RegTitle, RegN)
print*, "  > Files: A=temp(ann mean,C), B=prec(ann total,mm)"
call GetParallelFiles
call GetOutFile

do XFile = 1, FileN
  call GetMidFileIndex
  
  call LoadGloVec (InFile (XFile),    MapIDLReg, RegVector, Silent=1)	! load 'in' file
  call LoadGloVec (MidFile(XMidFile), MapIDLReg, MidVector, Silent=1)  	! load 'mid' file
  
  allocate (OutVector(RegN),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertBiome: Allocation failure: OV #####"
  OutVector=0	! initialised to 'undefined'
  BiomeReg=0	! number of allocations to each biome
  BiomeOther=0 ; BiomeOut=0
  
  do XReg = 1, RegN
    if (RegVector(XReg).NE.MissVal.AND.MidVector(XReg).NE.MissVal) then
      if (RegVector(XReg).GE.  30) RegVector(XReg)=  29.5
      if (RegVector(XReg).LT. -15) RegVector(XReg)= -14.5
      if (MidVector(XReg).GE.4500) MidVector(XReg)=4450
	
      XTmp=floor(RegVector(XReg))+16
      XPre=floor(MidVector(XReg)/100.0)+1
      OutVector(XReg)=Biome(XTmp,XPre)
      BiomeReg(Biome(XTmp,XPre))=BiomeReg(Biome(XTmp,XPre))+1
    else
      OutVector(XReg)=MissVal
    end if
  end do
  
  print "(a,16i6)", "   > ", (BiomeReg(XTmp),XTmp=1,16)
  write (99,"(i4,8i10)"), XFile, (BiomeReg(XTmp),XTmp=1,8)
  write (99,"(i4,8i10)"), XFile, (BiomeReg(XTmp),XTmp=9,16)
  
  SaveTitle = 'biomes'
  call SaveGlo (GridLongN,GridLatN,RegN,GridFile,OutFile(XFile),SaveTitle,OutVector,MapIDLReg)
  
  deallocate (RegVector,MidVector,OutVector, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertBiome: Deallocation failure: RV #####"
end do  

deallocate (MapIDLReg,RegSizes,RegNames,InFile,OutFile,MidFile,Biome,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ConvertBiome: Deallocation failure: file #####"

end subroutine ConvertBiome

!*******************************************************************************
! 13. report grid-box

subroutine ReportGridBox

print*, "  > Specify the common grid. "
call GridSelect (GridChosen,GridTitle,GridLongN,GridLatN,GridDataN,GridFile)
call RegSelect (GridChosen,GridLongN,GridLatN,GridDataN,MapIDLReg, RegSizes, RegNames, RegTitle, RegN)

print*, "  > Select the box reference (exe,wye): "
do
      read (*,*,iostat=ReadStatus), XExe,XWye
      if (XExe.LT.1.OR.XWye.LT.1.OR.XExe.GT.GridLongN.OR.XWye.GT.GridLatN) then
        ReadStatus=1 ; print*, "  > Out-of-range. Retry. "
      else if (MapIDLReg(XExe,XWye).EQ.MissVal) then
        ReadStatus=1 ; print*, "  > Not part of a region. Retry. "
      end if
      if (ReadStatus.LE.0) exit
end do

print*, "  > Select the files."
call GetBatch (Blank,InFile)
FileN  = size (InFile, 1)

do XFile=1,FileN
  GivenFile="" ; GivenFile = InFile(XFile)
  LastSlash = index(GivenFile,"/",.TRUE.)
  NameOnly  = adjustl(trim(GivenFile((LastSlash+1):80)))
  
  call LoadGloVec (InFile(XFile), MapIDLReg, RegVector, Silent=1)
  
  print "(a,i4,2a,f10.2)", "   > Batch item: ", XFile, ": ", trim(NameOnly), &
  		RegVector(MapIDLReg(XExe,XWye))
  write (99,"(a,i4,2a,f10.2)"), "   > Batch item: ", XFile, ": ", trim(NameOnly), &
  		RegVector(MapIDLReg(XExe,XWye))
end do

end subroutine ReportGridBox

!*******************************************************************************
! get 'mid' file index, depending on replication

subroutine GetMidFileIndex

if      (QReplicate.EQ.0) then		! work out 'mid' file index
    XMidFile = XFile
else if (QReplicate.EQ.1) then
    XMidFile = nint(mod(real(XFile),real(MidFileN)))
    if (XMidFile.EQ.0) XMidFile = MidFileN
else if (QReplicate.EQ.2) then
    if (XFile.EQ.1) XMidFile=0
    if (mod(real(XFile),(real(FileN)/real(MidFileN))).EQ.1) XMidFile=XMidFile+1
end if

if (QReplicate.GT.0) then
    print*
    GivenFile = InFile(XFile)
    LastSlash = index(GivenFile,"/",.TRUE.)
    NameOnly  = adjustl(trim(GivenFile((LastSlash+1):80)))
    print "(a,i4,2a)", "   > File A: ", XFile, ": ", trim(NameOnly)
    GivenFile = MidFile(XMidFile)
    LastSlash = index(GivenFile,"/",.TRUE.)
    NameOnly  = adjustl(trim(GivenFile((LastSlash+1):80)))
    print "(a,i4,2a)", "   > File B: ", XFile, ": ", trim(NameOnly)
else
    GivenFile = OutFile(XFile)
    LastSlash = index(GivenFile,"/",.TRUE.)
    NameOnly  = adjustl(trim(GivenFile((LastSlash+1):80)))
    print "(a,i4,2a)", "   > Batch item: ", XFile, ": ", trim(NameOnly)
end if
  
end subroutine GetMidFileIndex

!*******************************************************************************
! get parallel files

subroutine GetParallelFiles

print*, "  > Select OLD files 'A'. There must be a common substring."
call GetBatch (Blank,InFile)
FileN  = size (InFile, 1)

print*, "  > Select OLD files 'B'. There must be a common substring."
do
  call GetBatch (Blank,MidFile)
  MidFileN = size(MidFile, 1)
  if (FileN.NE.MidFileN) then
    if (mod(real(FileN),real(MidFileN)).EQ.0) then
      print*, "  > Replicate the files to match A ? (0=no,1,2=yes)"
      print*, "  >   1=(a:a,n+a,...) 2=(a:n*a+1,n*a+2,...)"
      do
        read (*,*,iostat=ReadStatus), QReplicate
        if (ReadStatus.LE.0.AND.QReplicate.GE.0.AND.QReplicate.LE.2) exit
      end do
    end if

    if (QReplicate.EQ.0) then
      print*, "  > Incorrect number of files. Try again."
      deallocate (MidFile, stat=AllocStat)
      if (AllocStat.NE.0) print*, "  > ##### ERROR: OpMultiGlo: Deallocation failure: MF #####"
    end if
  end if
  if (FileN.EQ.MidFileN.OR.QReplicate.GT.0) exit
end do

end subroutine GetParallelFiles

!*******************************************************************************
! get output file vector

subroutine GetOutFile

print*, "  > Enter the substring in the OLD A files to alter in the NEW files: "
do
	read (*,*,iostat=ReadStatus), InSub
	if (ReadStatus.GT.0) then	
			print*, "  > Bad format. Try again."
	else if (InSub.EQ."") then
			print*, "  > Blank not permitted. Try again."
	end if
	if (ReadStatus.LE.0.AND.InSub.NE."") exit
end do
InSub=trim(adjustl(InSub))

print "(3a)", "   > Enter the substring in the NEW files to replace '", trim(InSub), "':"
do
	read (*,*,iostat=ReadStatus), OutSub
	if (ReadStatus.GT.0) then	
			print*, "  > Bad format. Try again."
	else if (OutSub.EQ."") then
			print*, "  > Blank not permitted. Try again."
	end if
	if (ReadStatus.LE.0.AND.OutSub.NE."") exit
end do

InSubLen = len(trim(InSub))
  
allocate (OutFile (FileN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GloBulk: OutFile: Allocation failure #####"
OutFile = ""

do XFile = 1, FileN
      GivenFile = InFile(XFile)
      InSubBeg  = index(GivenFile,InSub(1:InSubLen))
      InFileLen = len(adjustl(trim(GivenFile)))
      
      OutFile(XFile) = GivenFile(1:(InSubBeg-1)) // trim(adjustl(OutSub)) // &
      		GivenFile((InSubBeg+InSubLen):InFileLen)
end do

end subroutine GetOutFile

!*******************************************************************************
! check output magnitudes

subroutine CheckOutMag

print*, "  > Restrict output magnitudes (0=no,1=yes) ?"
do
      read (*,*,iostat=ReadStatus), QOutMag
      if (ReadStatus.LE.0.AND.QOutMag.GE.0.AND.QOutMag.LE.1) exit
end do
if (QOutMag.EQ.1) then
  print*, "  > Magnitudes smaller than X are set to X. Enter X."
  do
      read (*,*,iostat=ReadStatus), MagMin
      if (ReadStatus.LE.0.AND.MagMin.GE.0) exit
  end do
  
  print*, "  > Magnitudes larger  than X are set to X. Enter X."
  do
      read (*,*,iostat=ReadStatus), MagMax
      if (ReadStatus.LE.0.AND.MagMax.GE.0) exit
  end do  
end if

end subroutine CheckOutMag

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

end program GloBulk