! initialmod.f90
! module procedure written by Tim Mitchell in Dec 1999
! last modification on 19.09.01
! includes all the routines necessary to start a job going
!	SetMissAccept, SeasonSelect, PeriodSelect, GridSelect, RawSelect, RegSelect, LoadRefReg, GetGrid

module InitialMod

use FileNames

implicit none

contains

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

subroutine SetMissAccept (MissAccept)

real, intent (out) :: MissAccept

integer :: ReadStatus

print*, "  > Enter the acceptable % of MissVal: "
do
	read (*,*,iostat=ReadStatus), MissAccept
	if (ReadStatus.LE.0 .AND. MissAccept.GT.0) exit
end do
	
end subroutine SetMissAccept

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

subroutine SeasonSelect (Month0,Month1,MonthN)

integer, intent (out) :: Month0,Month1,MonthN

integer :: GetMonth0, GetMonth1, GetMonthN
integer :: ReadStatus

print*, "  > Enter the first and last month nos. to analyse (overlaps OK, total>1): "
do
	read (*,*,iostat=ReadStatus), GetMonth0, GetMonth1
	if (ReadStatus.LE.0 .AND. GetMonth0.GE.1 .AND. GetMonth0.LE.12 &
				.AND. GetMonth1.GE.1 .AND. GetMonth1.LE.12 &
				.AND. GetMonth1.NE.GetMonth0) exit
end do

if ((GetMonth1-GetMonth0).GE.0) then
  GetMonthN = GetMonth1 - GetMonth0 + 1
else
  GetMonthN = 13 - GetMonth0 + GetMonth1
end if

Month0 = GetMonth0
Month1 = GetMonth1
MonthN = GetMonthN

end subroutine SeasonSelect

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

subroutine PeriodSelect (JobYearN, JobDecN, JobADYear)

integer, dimension (:), pointer :: JobADYear

integer, intent (out) 	:: JobYearN, JobDecN

real, parameter :: MissVal = -999.0

integer :: JobADYear0, JobADYear1
integer :: ADYear0, ADYear1
integer :: ReadStatus
integer :: AllocStat
integer :: XADYear, ArrayYear

print*, "  > Enter the first and last years AD to analyse: "
do
	read (*,*,iostat=ReadStatus), ADYear0, ADYear1
	if (ReadStatus.LE.0 .AND. ADYear1.GT.ADYear0) exit
end do
	
JobADYear0 =  10 * floor   (real(ADYear0) / 10.0)
JobADYear1 = (10 * ceiling (real(ADYear1) / 10.0)) - 1

JobYearN = JobADYear1 - JobADYear0 + 1
JobDecN  = JobYearN / 10

allocate (JobADYear (JobYearN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PeriodSelect: Allocation failure #####"

if (JobADYear0.LT.ADYear0) then
  do XADYear = JobADYear0, (ADYear0 - 1)
    ArrayYear = XADYear - JobADYear0 + 1
    JobADYear (ArrayYear) = MissVal
  end do
end if

do XADYear = ADYear0, ADYear1
    ArrayYear = XADYear - JobADYear0 + 1
    JobADYear (ArrayYear) = XADYear
end do

if (JobADYear1.GT.ADYear1) then
  do XADYear = (ADYear1 + 1), JobADYear1
    ArrayYear = XADYear - JobADYear0 + 1
    JobADYear (ArrayYear) = MissVal
  end do
end if

end subroutine PeriodSelect

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

subroutine GetDims (Grid,LongN,LatN)

integer, intent(in) :: Grid
integer, intent(out):: LongN, LatN
integer :: ReadStatus

if      (Grid.EQ. 2) then
  LongN = 128 ; LatN =  64
else if (Grid.EQ. 3) then
  LongN =  96 ; LatN =  73
else if (Grid.EQ. 5) then
  LongN =  96 ; LatN =  48
else if (Grid.EQ. 6) then
  LongN =  64 ; LatN =  56
else if (Grid.EQ.12) then
  LongN = 720 ; LatN = 360
else if (Grid.EQ.13) then
  LongN = 480 ; LatN = 313
else if (Grid.EQ.14) then
  LongN = 160 ; LatN = 106
else if (Grid.EQ.15) then
  LongN = 258 ; LatN = 228
else if (Grid.EQ.17) then
  LongN =  80 ; LatN =  72
else if (Grid.EQ.18) then
  LongN = 144 ; LatN =  73
else if (Grid.EQ.19) then
  LongN =  97 ; LatN =  73
else if (Grid.EQ.20) then
  LongN = 180 ; LatN = 290
else if (Grid.EQ.21) then
  LongN = 130 ; LatN = 242
else if (Grid.EQ.22) then
  LongN =  72 ; LatN =  36
else
  print*, "  > GetDims: your chosen grid is not specified."
  print*, "  > Enter the number of longs, lats on the grid: "
  do
	read (*,*,iostat=ReadStatus), LongN, LatN
	if (ReadStatus.LE.0 .AND. LongN.GT.0 .AND. LatN.GT.0) exit
  end do
end if

end subroutine GetDims

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

function GetGrid (LongN,LatN)

integer :: LongN,LatN,GetGrid			
integer	:: GridChosen,GridLongN,GridLatN,GridDataN			
character (len=10) :: GridTitle			
character (len=80) :: GridFile			

if      (LongN.EQ. 96.AND.LatN.EQ. 48) then
  GetGrid = 5
else if (LongN.EQ.128.AND.LatN.EQ. 64) then
  GetGrid = 2
else if (LongN.EQ. 64.AND.LatN.EQ. 56) then
  GetGrid = 6
else if (LongN.EQ.720.AND.LatN.EQ.360) then
  GetGrid = 12
else if (LongN.EQ.480.AND.LatN.EQ.313) then
  GetGrid = 13
else if (LongN.EQ.160.AND.LatN.EQ.106) then
  GetGrid = 14
else if (LongN.EQ.258.AND.LatN.EQ.228) then
  GetGrid = 15
else if (LongN.EQ. 80.AND.LatN.EQ. 72) then
  GetGrid = 17
else if (LongN.EQ.144.AND.LatN.EQ. 73) then
  GetGrid = 18
else if (LongN.EQ. 97.AND.LatN.EQ. 73) then
  GetGrid = 19
else if (LongN.EQ.180.AND.LatN.EQ.290) then
  GetGrid = 20
else if (LongN.EQ.130.AND.LatN.EQ.242) then
  GetGrid = 21
else if (LongN.EQ. 72.AND.LatN.EQ. 36) then
  GetGrid = 22
else
  call GridSelect (GridChosen,GridTitle,GridLongN,GridLatN,GridDataN,GridFile)
  GetGrid = GridChosen
end if

end function GetGrid

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

subroutine GridSelect (GridChosen,GridTitle,GridLongN,GridLatN,GridDataN,GridFile,CallGrid)

integer, intent (in), optional		:: CallGrid		! get grid details if grid already chosen		
integer, intent (out)			:: GridChosen,GridLongN,GridLatN,GridDataN			
character (len=10), intent (out) 	:: GridTitle			
character (len=80), intent (out) 	:: GridFile			

real, parameter :: MissVal = -999.0

integer :: ReadStatus				! status of user input request

if (present(CallGrid)) then
  GridChosen = CallGrid
else
  print*, "  > Select a grid (0=list): "	
  do
	read (*,*,iostat=ReadStatus), GridChosen
	if (ReadStatus.LE.0 .AND. GridChosen.EQ.0) then
		print*, "  >  1. obs (Hulme)	 2. ECHam4-OPYC3"
		print*, "  >  3. HadCM2		 4. HadCM3"
		print*, "  >  5. CGCM1           6. CSIRO"
		print*, "  >  7. CCSR-NIES	 8. ECHam3-LSG"
		print*, "  >  9. GFDL		10. NCAR-CGCM"
		print*, "  > 11. NCAR-PCM3	12. 0.5 deg"
		print*, "  > 13. 10' Europe	14. 0.5 Europe"
		print*, "  > 15. 10' ATEAM      16. **ReGrid**"
		print*, "  > 17. 10' Britain    18. Richard Jones"
		print*, "  > 19. h2grid(regs=3) 20. 5km NatGrid"	! any (19) regions are got from (3)
		print*, "  > 21. 5km NG cut     22. 5.0 deg"
        end if
	if (ReadStatus.LE.0 .AND. GridChosen.GE.1 .AND. GridChosen.LE.22) exit
  end do
end if

if (GridChosen.EQ.1)  GridFile = "./../../../constants/goglo/obshulme.ref"
if (GridChosen.EQ.2)  GridFile = "./../../../constants/goglo/echam4-opyc3.ref"
if (GridChosen.EQ.3)  GridFile = "./../../../constants/goglo/hadcm2.ref"
if (GridChosen.EQ.4)  GridFile = "./../../../constants/goglo/hadcm3.ref"
if (GridChosen.EQ.5)  GridFile = "./../../../constants/goglo/cgcm1.ref"
if (GridChosen.EQ.6)  GridFile = "./../../../constants/goglo/csiro.ref"
if (GridChosen.EQ.7)  GridFile = "./../../../constants/goglo/ccsr-nies.ref"
if (GridChosen.EQ.8)  GridFile = "./../../../constants/goglo/echam3-lsg.ref"
if (GridChosen.EQ.9)  GridFile = "./../../../constants/goglo/gfdl.ref"
if (GridChosen.EQ.10) GridFile = "./../../../constants/goglo/ncar-cgcm.ref"
if (GridChosen.EQ.11) GridFile = "./../../../constants/goglo/ncar-pcm3.ref"
if (GridChosen.EQ.12) GridFile = "./../../../constants/goglo/half-degree.ref"
if (GridChosen.EQ.13) GridFile = "./../../../constants/goglo/tenmin-euro.ref"
if (GridChosen.EQ.14) GridFile = "./../../../constants/goglo/half-degree-euro.ref"
if (GridChosen.EQ.15) GridFile = "./../../../constants/goglo/tenmin-ateam.ref"
if (GridChosen.EQ.17) GridFile = "./../../../constants/goglo/tenmin-brit.ref"
if (GridChosen.EQ.18) GridFile = "./../../../constants/goglo/richjones.ref"
if (GridChosen.EQ.19) GridFile = "./../../../constants/goglo/h2grid.ref"
if (GridChosen.EQ.20) GridFile = "./../../../constants/goglo/ng5.ref"
if (GridChosen.EQ.21) GridFile = "./../../../constants/goglo/ng5cut.ref"
if (GridChosen.EQ.22) GridFile = "./../../../constants/goglo/five-degree.ref"

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

read (4, fmt="(A10)"), GridTitle
read (4, fmt="(2(I4),I6)"), GridLongN, GridLatN, GridDataN

close (4)

end subroutine GridSelect

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

subroutine RawSelect (GridChosen,GridLongN,GridLatN,MapIDLReg,MapIDLRaw,MapRawReg)

integer, pointer, dimension (:,:) 	:: MapIDLReg,MapIDLRaw	
						! IDL mapping arrays, 1st=given,2nd=returned
integer, pointer, dimension (:) 	:: MapRawReg	
						! IDL mapping arrays, returned
integer, intent (in)			:: GridChosen, GridLongN, GridLatN			
						! Grid chosen
real, parameter :: MissVal = -999.0

integer :: AllocStat, ReadStatus
integer :: XLong, XLat, XEight
integer :: EightN, Long0, Long1
integer :: RawChosen, GridDataN

character (len=80) :: RawFile

if      (GridChosen.EQ. 3) then
  RawFile = "./../../../constants/goglo/hadcm2-raw.ref"
else if (GridChosen.EQ. 4) then
  RawFile = "./../../../constants/goglo/hadcm3-raw.ref"
else if (GridChosen.EQ.12) then
  print*, "  > Select the structure in which the raw data is stored:"
  print*, "  >  1. Link/DDC : N->S Gr->E"
  print*, "  >  2. IDL      : S->N DL->E"

  do
	read (*,*,iostat=ReadStatus), RawChosen
	if (ReadStatus.LE.0 .AND. RawChosen.GE.1 .AND. RawChosen.LE.2) exit
  end do
  
  if (RawChosen.EQ.1) RawFile = "./../../../constants/goglo/half-degree-raw-ddc.ref"
  if (RawChosen.EQ.2) RawFile = "./../../../constants/goglo/half-degree-raw-idl.ref"
else
  print*, "  > ##### RawSelect has no raw .ref for this grid. #####"
end if

GridDataN = GridLongN * GridLatN

allocate (MapIDLRaw (GridLongN,GridLatN), &
	  MapRawReg (GridDataN),          stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: RawSelect: Allocation failure #####"
MapIDLRaw = MissVal
MapRawReg = MissVal

EightN = GridLongN / 8

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

do XLat = 1, GridLatN
  do XEight = 1, EightN
    Long0 = ((XEight - 1) * 8) + 1
    Long1 =  Long0 + 7
    
    read (4, fmt="(8(I6))"), (MapIDLRaw (XLong,XLat), XLong=Long0,Long1)
  end do
end do

close (4)

do XLong = 1, GridLongN
  do XLat = 1, GridLatN
    if (MapIDLReg (XLong,XLat).NE.MissVal) then
      MapRawReg (MapIDLRaw (XLong,XLat)) = MapIDLReg (XLong, XLat)
    end if
  end do
end do

end subroutine RawSelect

!*******************************************************************************
! modified to allow loading of H2 regions onto h2grid (grid=19)
! MapIDLRaw and MapRawReg transferred to SelectRaw on 4.10.00

subroutine RegSelect (GridChosen, LongN, LatN, DataN, MapIDLReg, RegSizes, RegNames, &
			RegTitle, RegN, RegFileIn, RegFileOut)

integer, pointer, dimension (:,:) 		:: MapIDLReg	! shape of IDL mapping arrays
integer, pointer, dimension (:) 		:: RegSizes	! region sizes + reg-->raw

character (len=20), pointer, dimension (:) 	:: RegNames	! names of individual regions

integer, intent (in) 	:: GridChosen				! if RegFile is set, GridChosen can be MissVal
integer, intent (in)	:: LongN, LatN, DataN			! these MUST be set

integer, intent (out)			:: RegN			! the number of regions
character (len=80),intent(in),optional 	:: RegFileIn		! the file of the region set
character (len=80),intent(out),optional :: RegFileOut		! the file of the region set
character (len=40), intent (out) 	:: RegTitle		! the name of the region set

character (len=80) :: OutRefFile, FilePath

real, parameter :: MissVal = -999.0

integer :: ExeN,WyeN

if (present(RegFileIn)) then
	FilePath = RegFileIn
	ExeN = LongN ; WyeN = LatN
else
	call RegSelectFile (GridChosen,FilePath)	! if Grid=19, RSF gets H2 regions
	call GetDims (GridChosen,ExeN,WyeN)	
	if (LongN.NE.ExeN.OR.LatN.NE.WyeN) then
    		ExeN = MissVal ; WyeN = MissVal
    		print*, "  > ##### ERROR: RegSelect: found and called dims differ ##### "
	end if
end if

call LoadRefReg (FilePath,MapIDLReg,RegSizes,RegNames,&
			ExeN=ExeN,WyeN=WyeN,OutRefFile=OutRefFile,OutRefTitle=RegTitle)

RegN = size(RegSizes,1)
if (present(RegFileOut)) RegFileOut = OutRefFile

end subroutine RegSelect

!*******************************************************************************
! allows the operator to manually select the region set file

subroutine RegSelectFile (GridChosen,FilePath)

integer, intent (in) 		:: GridChosen 		! Grid chosen
integer :: RegionChosen				! internal use only
integer :: ReadStatus				! status of user input request

character (len=80), intent(out) :: FilePath	! file path of region details 
character (len=80) :: GivenFile			! file path of region details 

print*, "  > Select a region set (0=list): "
	
if      (GridChosen.EQ.2) then  
  do
	read (*,*,iostat=ReadStatus), RegionChosen
	
	if (RegionChosen.EQ.0) then
  		print*, "  >   4. One region per grid box"
	end if
	
	if (ReadStatus.LE.0 .AND. RegionChosen.GE.4 .AND. RegionChosen.LE.4) exit
  end do
  
  if (RegionChosen.EQ.4) FilePath = "./../../../constants/goglo/echam4-rpb.ref"
else if (GridChosen.EQ.3.OR.GridChosen.EQ.19) then  
  do
	read (*,*,iostat=ReadStatus), RegionChosen
	if (ReadStatus.LE.0 .AND. RegionChosen.EQ.0 ) then
    		print*, "  >   1. Land-Sea mask        'lsm'"
    		print*, "  >   2. North v South        'nvs'"
    		print*, "  >   3. Single global region 'sgg'"
    		print*, "  >   4. Region per box       'rpb'"
    		print*, "  >   5. Region per land box  'lan'"
    		print*, "  >   6. Countries            'cty'"
    		print*, "  >   7. IPCC regions         'ipc'"
    		print*, "  >   8. Topographic 25 box   't25'"
    		print*, "  >   9. Giorgi 2000          'gio'"
    		print*, "  >  10. Geometric 25 box     'g25'"
    		print*, "  >  11. Geometric 81 box     'g81'"
    		print*, "  >  12. Climatic             'cli'"
    		print*, "  >  13. Antarctica rpb       'ant'"
    		print*, "  >  14. Rest of World rpb    'row'"
    		print*, "  > 100. User specified       'usr'"
	end if
	if (ReadStatus.LE.0 .AND. RegionChosen.GE.1 .AND. RegionChosen.LE.14) exit
	if (ReadStatus.LE.0 .AND. RegionChosen.EQ.100) exit
  end do
  
  if (RegionChosen.EQ.1)  FilePath = "./../../../constants/goglo/hadcm2-landsea.ref"
  if (RegionChosen.EQ.2)  FilePath = "./../../../constants/goglo/hadcm2-northsouth.ref"
  if (RegionChosen.EQ.3)  FilePath = "./../../../constants/goglo/hadcm2-singleglo.ref"
  if (RegionChosen.EQ.4)  FilePath = "./../../../constants/goglo/hadcm2-regperbox.ref"
  if (RegionChosen.EQ.5)  FilePath = "./../../../constants/goglo/hadcm2-lan.ref"  
  if (RegionChosen.EQ.6)  FilePath = "./../../../constants/goglo/hadcm2-countries.ref"  
  if (RegionChosen.EQ.7)  FilePath = "./../../../constants/goglo/hadcm2-ipcc.ref"  
  if (RegionChosen.EQ.8)  FilePath = "./../../../constants/goglo/hadcm2-topog-25.ref"  
  if (RegionChosen.EQ.9)  FilePath = "./../../../constants/goglo/hadcm2-giorgi.ref"  
  if (RegionChosen.EQ.10) FilePath = "./../../../constants/goglo/hadcm2-geo-25.ref"  
  if (RegionChosen.EQ.11) FilePath = "./../../../constants/goglo/hadcm2-geo-81.ref"  
  if (RegionChosen.EQ.12) FilePath = "./../../../constants/goglo/hadcm2-clim.ref"  
  if (RegionChosen.EQ.13) FilePath = "./../../../constants/goglo/hadcm2-ant.ref"  
  if (RegionChosen.EQ.14) FilePath = "./../../../constants/goglo/hadcm2-row.ref"  

  if (RegionChosen.EQ.100) then  
    print*, "  > Enter path/name for region set .ref file: "
    do
	do
		read (*,*,iostat=ReadStatus), GivenFile
		if (ReadStatus.LE.0) exit
	end do
	
	inquire (file=GivenFile, name=FilePath)
	open (1, file=FilePath, status="old", iostat=ReadStatus)
	if (ReadStatus .EQ. 0) close (1)
	if (ReadStatus .EQ. 0) exit
    end do
  end if
else if (GridChosen.EQ.4) then  
  do
	read (*,*,iostat=ReadStatus), RegionChosen
	
	if (RegionChosen.EQ.0) then
  		print*, "  >   3. Single global region"
  		print*, "  >   4. One region per grid box"
	end if
	
	if (ReadStatus.LE.0 .AND. RegionChosen.GE.3 .AND. RegionChosen.LE.4) exit
  end do
  
  if (RegionChosen.EQ.3) FilePath = "./../../../constants/goglo/hadcm3-singleglo.ref"
  if (RegionChosen.EQ.4) FilePath = "./../../../constants/goglo/hadcm3-regperbox.ref"
else if (GridChosen.EQ.5) then  
  do
	read (*,*,iostat=ReadStatus), RegionChosen
	
	if (RegionChosen.EQ.0) then
  		print*, "  >   3. Single global region"
  		print*, "  >   4. One region per grid box"
	end if
	
	if (ReadStatus.LE.0 .AND. RegionChosen.GE.3 .AND. RegionChosen.LE.4) exit
  end do
  
  if (RegionChosen.EQ.3) FilePath = "./../../../constants/goglo/cgcm1-sgg.ref"
  if (RegionChosen.EQ.4) FilePath = "./../../../constants/goglo/cgcm1-rpb.ref"
else if (GridChosen.EQ.12) then
  
  do
	read (*,*,iostat=ReadStatus), RegionChosen
	
	if (RegionChosen.EQ.0) then
  		print*, "  >   2. North v South        'nvs'"
    		print*, "  >   3. Single global region 'sgg'"
  		print*, "  >   4. One region per box   'rpb'"
 	 	print*, "  >   5. Region per land box  'lan'"
  		print*, "  >   6. Countries v2        'cty2'"
  		print*, "  >   7. Countries v3 all-OK 'cty3'"
  		print*, "  >   8. E African Highlands  'eah'"
  		print*, "  >   9. Isothurm cities v1   'iso'"
  		print*, "  >  10. Koeppen area 1961-90 'kop'"
  		print*, "  >  11. Continents           'cnt'"
	end if
	
	if (ReadStatus.LE.0 .AND. RegionChosen.GE.2 .AND. RegionChosen.LE.11) exit
  end do
  
  if (RegionChosen.EQ.2) FilePath = "./../../../constants/goglo/half-degree-northsouth.ref"
  if (RegionChosen.EQ.3) FilePath = "./../../../constants/goglo/half-degree-sgg.ref"
  if (RegionChosen.EQ.4) FilePath = "./../../../constants/goglo/half-degree-rpb.ref"
  if (RegionChosen.EQ.5) FilePath = "./../../../constants/goglo/half-degree-lan.ref"
  if (RegionChosen.EQ.6) FilePath = "./../../../constants/goglo/half-degree-cty-v2.ref"
  if (RegionChosen.EQ.7) FilePath = "./../../../constants/goglo/half-degree-cty-v3.ref"
  if (RegionChosen.EQ.8) FilePath = "./../../../constants/goglo/half-degree-eah.ref"
  if (RegionChosen.EQ.9) FilePath = "./../../../constants/goglo/half-degree-isothurm-v1.ref"
  if (RegionChosen.EQ.10) FilePath = "./../../../constants/goglo/half-degree-kop.ref"
  if (RegionChosen.EQ.11) FilePath = "./../../../constants/goglo/half-degree-cnt.ref"
else if (GridChosen.EQ.13) then
  print*, "  >   1. Land-Sea mask        'lsm'"
  print*, "  >   4. One region per box   'rpb'"
  print*, "  >   5. Region per land box  'lan'"
  print*, "  > 100. User specified       'usr'"
  
  do
	read (*,*,iostat=ReadStatus), RegionChosen
	
	if (RegionChosen.EQ.0) then
  		print*, "  >   1. Land-Sea mask        'lsm'"
  		print*, "  >   4. One region per box   'rpb'"
  		print*, "  >   5. Region per land box  'lan'"
  		print*, "  > 100. User specified       'usr'"
	end if
	
	if (ReadStatus.LE.0 .AND. RegionChosen.GE.4 .AND. RegionChosen.LE.5) exit
	if (ReadStatus.LE.0 .AND. RegionChosen.EQ.1) 			     exit
	if (ReadStatus.LE.0 .AND. RegionChosen.EQ.100) 			     exit
  end do
  
  if (RegionChosen.EQ.1) FilePath = "./../../../constants/goglo/tenmin-euro-lsm.ref"
  if (RegionChosen.EQ.4) FilePath = "./../../../constants/goglo/tenmin-euro-rpb.ref"
  if (RegionChosen.EQ.5) FilePath = "./../../../constants/goglo/tenmin-euro-lan.ref"

  if (RegionChosen.EQ.100) then  
    print*, "  > Enter path/name for region set .ref file: "
    do
	do
		read (*,*,iostat=ReadStatus), GivenFile
		if (ReadStatus.LE.0) exit
	end do
	
	inquire (file=GivenFile, name=FilePath)
	open (1, file=FilePath, status="old", iostat=ReadStatus)
	if (ReadStatus .EQ. 0) close (1)
	if (ReadStatus .EQ. 0) exit
    end do
  end if
else if (GridChosen.EQ.15) then
  do
	read (*,*,iostat=ReadStatus), RegionChosen
	
	if (RegionChosen.EQ.0) then
  		print*, "  >   3. All-in-one"  
  		print*, "  >   4. Region per box"  
  		print*, "  >   5. Region per land box"  
	end if
	
	if (ReadStatus.LE.0 .AND. RegionChosen.GE.3 .AND. RegionChosen.LE.5) exit
  end do
  
  if (RegionChosen.EQ.3) FilePath = "./../../../constants/goglo/tenmin-ateam-one.ref"
  if (RegionChosen.EQ.4) FilePath = "./../../../constants/goglo/tenmin-ateam-rpb.ref"
  if (RegionChosen.EQ.5) FilePath = "./../../../constants/goglo/tenmin-ateam-lan.ref"
else if (GridChosen.EQ.17) then
  do
	read (*,*,iostat=ReadStatus), RegionChosen
	
	if (RegionChosen.EQ.0) then
  		print*, "  >   5. One region per land box"  
	end if
	
	if (ReadStatus.LE.0 .AND. RegionChosen.GE.5 .AND. RegionChosen.LE.5) exit
  end do
  
  if (RegionChosen.EQ.5) FilePath = "./../../../constants/goglo/tenmin-brit-lan.ref"
else if (GridChosen.EQ.18) then  
  do
	read (*,*,iostat=ReadStatus), RegionChosen
	
	if (RegionChosen.EQ.0) then
  		print*, "  >   4. One region per box"
  		print*, "  >  15. UK only"
  		print*, "  >  16. West Africa (waf)"
	end if
	
	if (ReadStatus.LE.0 .AND. RegionChosen.GE.4 .AND. RegionChosen.LE.16) exit
  end do
  
  if (RegionChosen.EQ. 4) FilePath = "./../../../constants/goglo/richjones-rpb.ref"
  if (RegionChosen.EQ.15) FilePath = "./../../../constants/goglo/richjones-uk.ref"
  if (RegionChosen.EQ.16) FilePath = "./../../../constants/goglo/richjones-waf.ref"
else if (GridChosen.EQ.20) then  
  do
	read (*,*,iostat=ReadStatus), RegionChosen
	
	if (RegionChosen.EQ.0) then
  		print*, "  >   4. One region per box"
	end if
	
	if (ReadStatus.LE.0 .AND. RegionChosen.GE.4 .AND. RegionChosen.LE.4) exit
  end do
  
  if (RegionChosen.EQ.4) FilePath = "./../../../constants/goglo/ng5-rpb.ref"
else if (GridChosen.EQ.21) then
  do
	read (*,*,iostat=ReadStatus), RegionChosen
	
	if (RegionChosen.EQ.0) then
  		print*, "  >   4. One region per box"
	end if
	
	if (ReadStatus.LE.0 .AND. RegionChosen.GE.4 .AND. RegionChosen.LE.4) exit
  end do
  
  if (RegionChosen.EQ.4) FilePath = "./../../../constants/goglo/ng5cut-rpb.ref"
else if (GridChosen.EQ.22) then
  do
	read (*,*,iostat=ReadStatus), RegionChosen
	
	if (RegionChosen.EQ.0) then
  		print*, "  >   4. One region per box"
	end if
	
	if (ReadStatus.LE.0 .AND. RegionChosen.GE.4 .AND. RegionChosen.LE.4) exit
  end do
  
  if (RegionChosen.EQ.4) FilePath = "./../../../constants/goglo/five-degree-rpb.ref"
else
  print*, "  > ### No regions listed in RegSelect for Grid ###"
end if

end subroutine RegSelectFile

!*******************************************************************************
! load .ref file describing region set
! takes a filepath, returns an array and two vectors

subroutine LoadRefReg (CallRefFile,MapIDLReg,RegSizes,RegNames,ExeN,WyeN,OutRefFile,OutRefTitle)

integer, pointer, dimension (:,:) 			:: MapIDLReg	! shape of IDL mapping arrays
integer, pointer, dimension (:) 			:: RegSizes	! region sizes + reg-->raw
character (len=20), pointer, dimension (:) 		:: RegNames	! names of individual regions

integer, intent(in), optional :: ExeN, WyeN	! optional: force MapIDLReg to ExeN,WyeN: RECOMMENDED ######
						! necessary when loading MapIDLReg for 97,73

character (len=80),intent (in) 		:: CallRefFile			! may be blank
character (len=80),intent(out),optional :: OutRefFile
character (len=40),intent(out),optional :: OutRefTitle
character (len=80) 			:: RefFile			! file path of region details 

real, parameter :: MissVal = -999.0

integer :: ReadStatus, AllocStat
integer :: RegN, LongEightN, LongTenN, LongN, LatN
integer :: XReg, XLat, XLong, XEight, XTen
integer :: Lat0, Lat1, Long0, Long1 
integer :: BlankTot, QContinue

character (len=40) :: RegTitle, LoadFormat

!***************************************
! load .ref describing region set

RefFile = LoadPath (CallRefFile,".ref")
if (present(OutRefFile)) OutRefFile = RefFile

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

read (4, fmt="(A40)"), 		RegTitle
read (4, fmt="(I6,2(I4))"),	RegN, LongN, LatN

if (present(OutRefTitle)) OutRefTitle = RegTitle

if (present(ExeN).AND.present(WyeN)) then
  if (LongN.NE.ExeN.OR.LatN.NE.WyeN) then
    if (ExeN.EQ.97.AND.LongN.EQ.96) then		! ensures that a 97,73 MapIDLReg is returned
      QContinue = 2
    else
      print*, "  > ##### ERROR: LoadRefReg: called and file dims do not match #####"
      QContinue = 0
    end if
  else
    QContinue = 1
  end if
else
  QContinue = 1
end if
  
if      (QContinue.EQ.1) then
  allocate (MapIDLReg (LongN, LatN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: RegSelect: Allocation failure: MapIDLReg 1 #####"
  MapIDLReg = MissVal
else if (QContinue.EQ.2) then
  allocate (MapIDLReg (ExeN, WyeN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: RegSelect: Allocation failure: MapIDLReg 2 #####"
  MapIDLReg = MissVal
end if

allocate (RegSizes  (RegN), &
	  RegNames  (RegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: RegSelect: Allocation failure:  Reg* #####"
RegSizes  = 0 ; RegNames  = ""

do XReg = 1, RegN
    read (4, fmt="(A20)"), RegNames (XReg)
end do

if ((real(LongN)/8.0).EQ.nint(real(LongN)/8.0)) then
  LoadFormat = "(8i6)"
else if ((real(LongN)/10.0).EQ.nint(real(LongN)/10.0)) then
  LoadFormat = "(10i6)"
else if ((real(LongN)/6.0).EQ.nint(real(LongN)/6.0)) then
  LoadFormat = "(6i6)"
else
  print*, "  > ##### ERROR: LoadRefReg: LongN not divisible by 10 or 8 or 6 #####"
end if

if      (trim(LoadFormat).EQ."(8i6)") then
 LongEightN = LongN / 8

 do XLat = 1, LatN
  do XEight = 1, LongEightN
    Long0 = ((XEight - 1) * 8) + 1
    Long1 =  Long0 + 7
    read (4, fmt=LoadFormat), (MapIDLReg (XLong,XLat), XLong=Long0,Long1)
  end do
 end do
else if (trim(LoadFormat).EQ."(10i6)") then
 LongTenN = LongN / 10

 do XLat = 1, LatN
  do XTen = 1, LongTenN
    Long0 = ((XTen - 1) * 10) + 1
    Long1 =  Long0 + 9
    read (4, fmt=LoadFormat), (MapIDLReg (XLong,XLat), XLong=Long0,Long1)
  end do
 end do
else if (trim(LoadFormat).EQ."(6i6)") then
 LongTenN = LongN / 6

 do XLat = 1, LatN
  do XTen = 1, LongTenN
    Long0 = ((XTen - 1) * 6) + 1
    Long1 =  Long0 + 5
    read (4, fmt=LoadFormat), (MapIDLReg (XLong,XLat), XLong=Long0,Long1)
  end do
 end do
end if

close (4)

!***************************************
! process region set

do XLong = 1, LongN
  do XLat = 1, LatN
    if (MapIDLReg (XLong,XLat).NE.MissVal) then
      RegSizes  (MapIDLReg (XLong,XLat)) = RegSizes (MapIDLReg (XLong,XLat)) + 1    
    end if
  end do
end do

BlankTot = 0
do XReg = 1, RegN
  if (RegSizes(XReg).EQ.0) then
    BlankTot = BlankTot + 1
  end if
end do
if (BlankTot.GT.0) print*, "  > ##### Blank regions found, totalling: ", BlankTot 

end subroutine LoadRefReg 

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

end module InitialMod