! makeregionref.f90
! main program written by Tim Mitchell on 07.12.99
! last modification on 15.12.99
! generates and saves a region .ref file
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
! 	-o ./../goglo/makeregionref filenames.f90 initialmod.f90 ./../goglo/makeregionref.f90

program MakeRegionRef

use FileNames
use InitialMod

implicit none

integer, dimension (:,:), allocatable 	:: MapIDLReg, LandSeaBin, NorthSouthBin

integer, pointer, dimension (:,:) :: MyMapIDLRaw

character (len=20), dimension (:), allocatable :: RegNames

character (len=80) :: FilePath, SaveFilePath, SelectFile, LandSeaFilePath, MyFilePath
character (len=10) :: MyName
character (len=10) :: LineFormat
character (len=40) :: Title

integer :: MyNo, MyLongN, MyLatN, MyDataN
integer :: RawLat, RawLong
integer :: RegN
integer	:: XLong, XLat, XDatum, XReg
integer :: Lat0, Lat1, LatStep
integer :: Long0, Long1
integer :: LongHalfTime0, LongHalfTime1
integer :: LongSelect, LatSelect
integer :: LatEightN, LongEightN, LongTenN, LongSixN, XEight, XTen, XSix
integer	:: AllocStat
integer :: MaskChoice
integer :: NorthSouth		! 1=northwards, 2=southwards
integer :: GreenDate		! start: 1=Greenwich, 2=DateLine
integer :: ReadStatus		! status of input from user
integer :: QType		! which type of ref file to generate
integer :: QNoPoles
integer :: HalfLat0, HalfLat1

real, parameter :: MissVal = -999.0

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

print*
print*, "  > MakeRegionRef: generate a model reference file"
print*

call GridSelect (MyNo, MyName, MyLongN, MyLatN, MyDataN, MyFilePath)

allocate (MapIDLReg (MyLongN, MyLatN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: allocation failure #####"
MapIDLReg = nint (MissVal)

print*, "  > Enter a choice: "
print*, "  >   1. Land-sea mask"
print*, "  >   2. N-S hemispheres"
print*, "  >   3. Single global region"
print*, "  >   4. One region per grid box"
print*, "  >   5. One region per land box only"
do
	read (*,*,iostat=ReadStatus), MaskChoice
	if (ReadStatus.LE.0) exit
end do
	
print*, "  > Exclude the first and last latitude bands ? (1=no,2=yes)"
do
	read (*,*,iostat=ReadStatus), QNoPoles
	if (ReadStatus.LE.0.AND.QNoPoles.GE.1.AND.QNoPoles.LE.2) exit
end do

if (MaskChoice.EQ.1) call SetLandSea
if (MaskChoice.EQ.2) call SetNorthSouth
if (MaskChoice.EQ.3) call SetSingleGlobe
if (MaskChoice.EQ.4) call SetRegPerBox
if (MaskChoice.EQ.5) call SetRegPerLandBox

call SaveRegionRef

print*

contains

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

subroutine GetLandSeaMask

allocate (LandSeaBin (MyLongN, MyLatN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: allocation failure #####"
LandSeaBin = -1

print*, "  > Enter the filepath of the DDC binary land-sea mask: " 
do
	do
		read (*,*,iostat=ReadStatus), SelectFile
		if (ReadStatus.LE.0) exit
	end do
	
	inquire (file=SelectFile, name=LandSeaFilePath)
	open (1, file=LandSeaFilePath, status="old", iostat=ReadStatus)
	if (ReadStatus .EQ. 0) close (1)
	if (ReadStatus .EQ. 0) exit
end do

print*, '  > Enter the line format e.g. "(96(I1))"'
do
	read (*,*,iostat=ReadStatus), LineFormat
	if (ReadStatus.LE.0) exit
end do

open (4, file=LandSeaFilePath, status="old", access="sequential", action="read")

do XLat = 1, MyLatN
   read (4, fmt=LineFormat), LandSeaBin (1:MyLongN,XLat)
end do

do XLong = 1, MyLongN
  do XLat = 1, MyLatN
    RawLong = MyMapIDLRaw (XLong, XLat)
    RawLat  = 1
    
    do 
      if (RawLong.LE.MyLongN) exit
      RawLong = RawLong - MyLongN
      RawLat  = RawLat  + 1
    end do
    
    MapIDLReg (XLong, XLat) = LandSeaBin (RawLong, RawLat) + 1
  end do
end do 

close (4)

deallocate (LandSeaBin)

end subroutine GetLandSeaMask

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

subroutine SetLandSea

RegN = 2
allocate (RegNames (RegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: allocation failure #####"
RegNames (1) = "Land"
RegNames (2) = "Sea"

call GetLandSeaMask

if (QNoPoles.EQ.2) then
  do XLong = 1, MyLongN
    MapIDLReg (XLong,1)      = nint (MissVal)
    MapIDLReg (XLong,MyLatN) = nint (MissVal)
  end do
end if

end subroutine SetLandSea

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

subroutine SetNorthSouth

RegN = 2

allocate (RegNames (RegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: allocation failure #####"

RegNames (1) = "North"
RegNames (2) = "South"

HalfLat0 = MyLatN / 2
HalfLat1 = HalfLat0 + 1

do XLat = 1, HalfLat0
  do XLong = 1, MyLongN
    MapIDLReg (XLong, XLat) = 2
  end do
end do

do XLat = HalfLat1, MyLatN
  do XLong = 1, MyLongN
    MapIDLReg (XLong, XLat) = 1
  end do
end do

if (QNoPoles.EQ.2) then
  do XLong = 1, MyLongN
    MapIDLReg (XLong,1)      = nint (MissVal)
    MapIDLReg (XLong,MyLatN) = nint (MissVal)
  end do
end if

end subroutine SetNorthSouth

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

subroutine SetSingleGlobe

RegN = 1
allocate (RegNames (RegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: allocation failure #####"
RegNames (1) = "Globe"

do XLong = 1, MyLongN
  do XLat = 1, MyLatN
    MapIDLReg (XLong, XLat) = 1
  end do
end do

if (QNoPoles.EQ.2) then
  do XLong = 1, MyLongN
    MapIDLReg (XLong,1)      = nint (MissVal)
    MapIDLReg (XLong,MyLatN) = nint (MissVal)
  end do
end if

print*, "  > Allocated single global region."

end subroutine SetSingleGlobe

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

subroutine SetRegPerBox

XReg = 0
do XLong = 1, MyLongN
  do XLat = 1, MyLatN
    if (QNoPoles.EQ.1) then
      XReg = XReg + 1
      MapIDLReg (XLong, XLat) = XReg
    else
     if (XLat.NE.1.AND.XLat.NE.MyLatN) then
      XReg = XReg + 1
      MapIDLReg (XLong, XLat) = XReg
     end if
    end if
  end do
end do
RegN = XReg

allocate (RegNames (RegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: allocation failure #####"
RegNames = "a-box"

end subroutine SetRegPerBox

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

subroutine SetRegPerLandBox

call GetLandSeaMask

XReg = 0
do XLong = 1, MyLongN
  do XLat = 1, MyLatN
   if (LandSeaBin (XLong, XLat) .EQ. 0) then
    if (QNoPoles.EQ.1) then
      XReg = XReg + 1
      MapIDLReg (XLong, XLat) = XReg
    else
     if (XLat.NE.1.AND.XLat.NE.MyLatN) then
      XReg = XReg + 1
      MapIDLReg (XLong, XLat) = XReg
     end if
    end if
   end if
  end do
end do
RegN = XReg

allocate (RegNames (RegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: allocation failure #####"
RegNames = "a-land-box"

end subroutine SetRegPerLandBox

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

subroutine SaveRegionRef

print*, "  > Enter the filepath of the region file: "
do
	do
		read (*,*,iostat=ReadStatus), SelectFile
		if (ReadStatus.LE.0) exit
	end do
	
	inquire (file=SelectFile, name=SaveFilePath)
	open (2, file=SaveFilePath, status="new", access="sequential", iostat=ReadStatus)
	if (ReadStatus .EQ. 0) close (2)
	if (ReadStatus .EQ. 0) exit
end do

print*, "  > Enter the title of the region file: "
do
	read (*,*,iostat=ReadStatus), Title
	if (ReadStatus.LE.0 .AND. Title.NE."") exit
end do
	
open (5, file=SaveFilePath, status="replace", access="sequential", action="write")

write (5, fmt="(A40)"), 	Title
write (5, fmt="(I6,2(I4))"),	RegN, MyLongN, MyLatN

do XReg = 1, RegN
    write (5, fmt="(A20)"), RegNames (XReg)
end do

if      (mod(real(MyLongN),8.0).EQ.0.0) then
 print*, "  > Saving in eights."
 LongEightN = MyLongN / 8
 do XLat = 1, MyLatN
  do XEight = 1, LongEightN
    Long0 = ((XEight - 1) * 8) + 1
    Long1 =  Long0 + 7
    write (5, fmt="(8(I6))"), MapIDLReg (Long0:Long1,XLat)
  end do
 end do
else if (mod(real(MyLongN),10.0).EQ.0.0) then
 print*, "  > Saving in tens."
 LongTenN = MyLongN / 10
 do XLat = 1, MyLatN
  do XTen = 1, LongTenN
    Long0 = ((XTen - 1) * 10) + 1
    Long1 =  Long0 + 9
    write (5, fmt="(10(I6))"), MapIDLReg (Long0:Long1,XLat)
  end do
 end do
else if (mod(real(MyLongN),6.0).EQ.0.0) then
 print*, "  > Saving in sixes."
 LongSixN = MyLongN / 6
 do XLat = 1, MyLatN
  do XSix = 1, LongSixN
    Long0 = ((XSix - 1) * 6) + 1
    Long1 =  Long0 + 5
    write (5, fmt="(6(I6))"), MapIDLReg (Long0:Long1,XLat)
  end do
 end do
else
  print*, "  > ##### ERROR: LongN not divisible by 10 or 8 or 6 #####"
end if

close (5)

end subroutine SaveRegionRef

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

end program MakeRegionRef