! makeborders.f90
! loads border data from arc file containing polys
! saves border data in .bnd form
! f90 -o makeborders filenames.f90 gridops.f90 makeborders.f90
! written by Tim Mitchell on 30.07.01

program MakeBorders

use FileNames
use GridOps

implicit none

integer, parameter :: nine = selected_real_kind (11,99)

real (kind=nine), dimension (4) :: Line

integer, allocatable, dimension (:,:)	:: PolyRef

real (kind=nine) :: Lat,Long,NodeX,NodeY,StartX,StartY

integer :: ReadStatus, AllocStat
integer :: HeadN,FootN,NodeN,PolyN,GapN,DoubleN,SumNodeN
integer :: XPoly,XHead,XFoot,XGap,XNode,XNodeLoad,XDatum
integer :: QTransform,QHeadFoot,QOrigForm

character (len=80), parameter :: Blank80 = ""
character (len=04), parameter :: Blank04 = "    "

character (len=80) :: LoadFile, PolyFile, HeadFile, DataFile, LoadFormat, SaveFile, Trash

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

open (99,file="/cru/u2/f709762/data/scratch/log-mkb.dat",status="replace",action="write")

print*
print*, "  > ***** MakeBorders: loads border data and saves as .bnd *****"
print*

call MakeChoices

if      (QOrigForm.EQ.1) then
  call LoadInfoPoly
  call LoadSavePoly
else if (QOrigForm.EQ.2) then
  call LoadSaveE00
end if

print*

close (99)

contains

!*******************************************************************************
! make choices

subroutine MakeChoices

print*, "  > Load a .poly (=1) or .e00 (=2) file ?"
do
	read (*,*,iostat=ReadStatus), QOrigForm
	if (ReadStatus.LE.0.AND.QOrigForm.GE.1.AND.QOrigForm.LE.2) exit
end do

if      (QOrigForm.EQ.1) then
  LoadFile = LoadPath (Blank80,"poly")
else if (QOrigForm.EQ.2) then
  LoadFile = LoadPath (Blank80,".e00")
end if

if      (QOrigForm.EQ.1) then
  LoadFormat = "(2f18.6)"
  HeadN = 1 ; FootN = 1
end if

print*, "  > Transform from National Grid to lat/long (1=no,2=yes) ?"
do
	read (*,*,iostat=ReadStatus), QTransform
	if (ReadStatus.LE.0.AND.QTransform.GE.1.AND.QTransform.LE.2) exit
end do

print*, "  > Specify the file to save."
SaveFile = SavePath (Blank80,".bnd")

end subroutine MakeChoices

!*******************************************************************************
! load .poly file info

subroutine LoadInfoPoly

print*, "  > Getting polygon information from load file..."

PolyFile = "/cru/u2/f709762/data/scratch/makeborders-polys.dat"

open  (1,file=LoadFile,status="old",action="read")
open  (2,file=PolyFile,status="scratch",action="readwrite")

NodeN = 0 ; PolyN = 0 ; GapN = 0 ; SumNodeN = 0
do
   read (1,fmt="(a80)",iostat=ReadStatus), Trash
   GapN = GapN + 1
   
   if      (adjustl(trim(Trash)).EQ."-99999") then
    do
       read (1,fmt="(a80)"), Trash
       GapN = GapN + 1
       if (adjustl(trim(Trash)).EQ."END") exit
    end do
   else if (adjustl(trim(Trash)).EQ."END") then
    ReadStatus = -1
   else if (ReadStatus.EQ.0) then
    if (HeadN.GT.1) then
      do XHead = 2, HeadN
        read (1,fmt="(a80)"), Trash
        GapN = GapN + 1
      end do
    end if
    
    read (1,fmt=LoadFormat), StartX, StartY
    NodeN = 1
    do
     read (1,fmt=LoadFormat), NodeX, NodeY
     NodeN = NodeN + 1
     if (NodeX.EQ.StartX.AND.NodeY.EQ.StartY) exit
    end do
    
    SumNodeN = SumNodeN + NodeN
    write ( 2,"(2i6)"), GapN, NodeN  	  
    write (99,"(3i10)"), GapN, NodeN, SumNodeN  	  
    GapN = 0 ; NodeN = 0 ; PolyN = PolyN + 1
    
    if (FootN.GT.0) then
      do XFoot = 1, FootN
        read (1,fmt="(a80)",iostat=ReadStatus), Trash
        GapN = GapN + 1
      end do
    end if
   end if
   
   if (ReadStatus.NE.0) exit
end do

close (1) ; rewind (2)

allocate (PolyRef (PolyN,2), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: MakeBorders: Allocation failure #####"

do XPoly = 1, PolyN
  read (2,"(2i6)"), PolyRef(XPoly,1), PolyRef(XPoly,2)	
end do

close (2)

end subroutine LoadInfoPoly

!*******************************************************************************
! load from original file and save to .bnd file

subroutine LoadSavePoly

print*, "  > Loading and saving data..."

open  (1,file=LoadFile,status="old",action="read")
open  (3,file=SaveFile,status="new",action="write")

write (3,"(a40,i6)"), "bounds file, containing total polygons: ", PolyN

do XPoly = 1, PolyN
  write (3,"(a13,2i6))"), "poly, nodes: ", XPoly, PolyRef(XPoly,2)
  
  do XHead = 1, PolyRef(XPoly,1)
    read  ( 1,fmt="(a80)"), Trash
  end do
  
  do XNode = 1, PolyRef(XPoly,2)
    read  ( 1,fmt=LoadFormat), NodeX,NodeY    
    
    if (QTransform.EQ.1) then
      Lat = NodeY ; Long = NodeX
    else
      call NatGridToLatLong (NodeX,NodeY,Lat,Long)
    end if
    
    write (3,"(2f12.4)"), Long,Lat
  end do
end do

close (3)
close (1)

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

end subroutine LoadSavePoly

!*******************************************************************************
! load from original file and save to .bnd file

subroutine LoadSaveE00

print*, "  > Loading and saving data..."

HeadFile = "/cru/u2/f709762/data/scratch/makeborders-head.dat"
DataFile = "/cru/u2/f709762/data/scratch/makeborders-data.dat"

open  (1,file=LoadFile,status="old",action="read")
open  (4,file=DataFile,status="replace",action="write")

read  (1,fmt="(a80)"), Trash			! get past headers
read  (1,fmt="(a80)"), Trash

XPoly = 0
do
  read  (1,fmt="(a60,i10)"), Trash, NodeN
  
  if (NodeN.GT.0) then
    XPoly = XPoly + 1

    write (4,fmt="(a13,2i6))"), "poly, nodes: ", XPoly, NodeN
  
    XNode = -1
    do
      XNode = XNode + 2

      if (XNode.EQ.NodeN) then
      	DoubleN = 1 ; LoadFormat = "(2e14.7)"
      else
      	DoubleN = 2 ; LoadFormat = "(4e14.7)"
      end if
      
      read (1,fmt=LoadFormat), (Line(XDatum),XDatum=1,(DoubleN*2))
      
      do XNodeLoad = 0, (DoubleN-1)
        NodeX = Line((XNodeLoad*2)+1) ; NodeY = Line((XNodeLoad*2)+2)
        
        if (QTransform.EQ.1) then
          Lat = NodeY ; Long = NodeX
        else
          call NatGridToLatLong (NodeX,NodeY,Lat,Long)
        end if
    
        write (4,"(2f12.4)"), Long,Lat
      end do
      
      if ((XNode+2).GT.NodeN) exit
    end do
  end if
  
  if (NodeN.EQ.0) exit
end do

close (4)
close (1)

open  (5,file=HeadFile,status="replace",action="write")
write (5,"(a40,i6)"), "bounds file, containing total polygons: ", XPoly
close (5)

call system ('cat ' // trim(HeadFile) // ' ' // trim(DataFile) // ' > ' // trim(SaveFile))
call system ('rm '  // trim(HeadFile))
call system ('rm '  // trim(DataFile))

end subroutine LoadSaveE00

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

end program MakeBorders