! gridops.f90
! written by Tim Mitchell, Tyndall Centre for Climate Change Research
! 27th July 2001, 14th September 2001
! module in which routines to handle transformations between grids are held
!	also routine to obtain distance between two points on globe
! contains:
!	NatGridToLatLong: converts from National Grid numbers to   lat/long
!	LatLongToNatGrid: converts to   National Grid numbers from lat/long
!	GetDistance	: gives distance between two points on globe

module GridOps

implicit none

contains

!*******************************************************************************
! supply with Got (NLat,NLon)=.TRUE. where a stn is held
! returns Need=.TRUE. for all boxes within range of >0 stns
! Got and Need must already be allocated ; Range is given in km

subroutine IdentifyRelevant (Got,Need,Range)

logical, pointer, dimension (:,:)	:: Got,Need
integer, allocatable, dimension (:)	:: ExeAdd
real, intent(in)			:: Range

real, parameter :: MissVal = -999.0

real    :: BoxLatDeg,BoxLonDeg, Extra,Dist, Lat,Lon
integer :: NLat,NLon,NExe,NWye, AllocStat
integer :: XLat,XLon,XExe,XWye, XSeekLat,XSeekLon

NLat=size(Got,1) ; NLon=size(Got,2)
BoxLatDeg=180.0/real(NLat) ; BoxLonDeg=360.0/real(NLon)
Need=.FALSE.

allocate (ExeAdd(NLat),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### Screen: dealloc #####"
ExeAdd=MissVal

Extra=0 ; NWye=MissVal		! find no. of N-S boxes to add
do
  Extra=Extra+1
  Dist=GetDistance(0.0,0.0,(Extra*BoxLatDeg),0.0)
  if (Dist.GT.Range) NWye=Extra-1
  if (NWye.EQ.0) NWye=1
  if (NWye.NE.MissVal.OR.(Extra*2).GE.NLat) exit
end do
if (NWye.EQ.MissVal) NWye=Extra

do XLat=1,NLat			! find no. of E-W boxes to add, by lat
  Extra=0 ; Lat=-180.0+((real(NLat)-0.5)*BoxLatDeg)
  do
    Extra=Extra+1
    Dist=GetDistance(Lat,0.0,Lat,(Extra*BoxLonDeg))
    if (Dist.GT.Range) ExeAdd(XLat)=Extra-1
    if (ExeAdd(XLat).EQ.0) ExeAdd(XLat)=1
    if (ExeAdd(XLat).NE.MissVal.OR.(Extra*2).GE.NLon) exit
  end do
  if (ExeAdd(XLat).EQ.MissVal) ExeAdd(XLat)=Extra
!  write (99,"(2i4,f9.2)"), XLat,ExeAdd(XLat),Lat	! @@@@@@@@@@@@@@@@
end do

do XLat=1,NLat
  NExe=ExeAdd(XLat)
  
  do XLon=1,NLon
    if (Got(XLat,XLon)) then
      do XExe=(0-NExe),NExe
        XSeekLon=XLon+XExe
	if (XSeekLon.LE.0) 	XSeekLon=XSeekLon+NLon
	if (XSeekLon.GT.NLon) 	XSeekLon=XSeekLon-NLon
	
	do XWye=(0-NWye),NWye
	  XSeekLat=XLat+XWye
	  if (XSeekLat.LE.0) 	XSeekLat=1-XSeekLat
	  if (XSeekLat.GT.NLat) XSeekLat=1-XSeekLat+(2*NLat)
	  
	  Need(XSeekLat,XSeekLon)=.TRUE.
	end do
      end do
    end if
  end do
end do

end subroutine IdentifyRelevant

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

function GetDistance (Lat1,Lon1,Lat2,Lon2)

real :: GetDistance
real :: Lat1,Lon1,Lat2,Lon2
real :: KayB,KayR
real :: VarThet1,VarThet2,VarPhi1,VarPhi2,VarX1,VarX2,VarY1,VarY2,VarZ1,VarZ2
real :: VarH,VarTemp,VarAngle

KayB=180.0/(22.0/7.0) ; KayR=6367.65

VarThet2=Lat2/KayB ; VarThet1=Lat1/KayB
VarPhi2=Lon2/KayB ; VarPhi1=Lon1/KayB

VarX1=sin(VarPhi1)*cos(VarThet1) ; VarX2=sin(VarPhi2)*cos(VarThet2)
VarY1=cos(VarPhi1)*cos(VarThet1) ; VarY2=cos(VarPhi2)*cos(VarThet2)
VarZ1=sin(VarThet1) ; VarZ2=sin(VarThet2)

VarH = (VarX1-VarX2)**2+(VarY1-VarY2)**2+(VarZ1-VarZ2)**2

VarTemp=(1+1-Varh)/2

if(VarTemp.lt.-1.0000000) VarTemp=-1.0000
if(VarTemp.gt. 1.0000000) VarTemp= 1.0000

VarAngle=acos(VarTemp)
GetDistance=VarAngle*KayR

end function GetDistance

!*******************************************************************************
! converts from National Grid numbers to lat/long
! East,North should be given in metres
! Lat,Long are returned in degrees as a decimal, rather than degrees/minutes/seconds

subroutine NatGridToLatLong (East,North,LatDeg,LonDeg)

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

real (kind=big), intent (in)  :: East,North
real (kind=big), intent (out) :: LatDeg,LonDeg

real (kind=big), parameter :: Ka = 6377563.396
real (kind=big), parameter :: Kb = 6356256.910
real (kind=big), parameter :: Kpi = 3.14159265358979
real (kind=big), parameter :: KN0 = -100000.0
real (kind=big), parameter :: KE0 =  400000.0
real (kind=big), parameter :: KF0 = 0.9996012717
real (kind=big), parameter :: KThi0 = 49.0
real (kind=big), parameter :: KLam0 = -2.0

real (kind=big) :: Lat,Lon
real (kind=big) :: Ve2, VThi0, VLam0, VThiDash, VM, Vro, Vv, Vnu2, Vsin, Vtan, Vsec
real (kind=big) :: V7, V8, V9, V10, V11, V12, V12A, Veas

integer :: QClose = 0

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

Ve2      = ((Ka*Ka)-(Kb*Kb))/(Ka*Ka)
VThi0    = (KThi0 / real(180,big)) * Kpi
VLam0    = (KLam0 / real(180,big)) * Kpi
VThiDash = ((North-KN0)/(Ka*KF0)) + VThi0

do
  VM = CalcM (Ka,Kb,KF0,VThi0,VThiDash)
  
  if ((North-KN0-VM).LT.real(0.001,big)) QClose = 1
  
  VThiDash = ((North-KN0-VM)/(Ka*KF0)) + VThiDash
  
  if (QClose.EQ.1) exit
end do

Vsin = sin(VThiDash) ; Vtan = tan(VThiDash) ; Vsec = real(1,big) / cos(VThiDash)

Vv   = Ka * KF0 * ((real(1,big)-(Ve2*Vsin*Vsin))**(real(-0.5,big)))

Vro  = Ka * KF0 * (real(1,big)-Ve2) * ((real(1,big)-(Ve2*Vsin*Vsin))**(real(-1.5,big)))

Vnu2 = (Vv/Vro) - real(1,big)

V7   = Vtan / (real(2,big)*Vro*Vv)

V8   = (Vtan / (real(24,big)*Vro*Vv*Vv*Vv)) * (real(5,big)+(real(3,big)*Vtan*Vtan)+Vnu2-(real(9,big)*Vtan*Vtan*Vnu2))

V9   = (Vtan / (real(720,big)*Vro*Vv*Vv*Vv*Vv*Vv)) * (real(61,big)+(real(90,big)*Vtan*Vtan)+(real(45,big)*Vtan*Vtan*Vtan*Vtan))

V10  = Vsec / Vv

V11  = (Vsec / (real(6,big)*Vv*Vv*Vv)) * ((Vv/Vro)+(real(2,big)*Vtan*Vtan))

V12  = (Vsec / (real(120,big)*Vv*Vv*Vv*Vv*Vv)) * (real(5,big)+(real(28,big)*Vtan*Vtan)+(real(24,big)*Vtan*Vtan*Vtan*Vtan))

V12A = (Vsec / (real(5040,big)*(Vv**7))) * (real(61,big)+(real(662,big)*(Vtan**2))+(real(1320,big)*(Vtan**4))+(real(720,big)*(Vtan**6)))

Veas = East - KE0

LatDeg = VThiDash - (V7*(Veas**2)) + (V8*(Veas**4)) - (V9*(Veas**6))

LonDeg = VLam0 + (V10*Veas) - (V11*(Veas**3)) + (V12*(Veas**5)) - (V12A*(Veas**7))

LatDeg = (LatDeg/Kpi) * real(180,big)

LonDeg = (LonDeg/Kpi) * real(180,big)

end subroutine NatGridToLatLong

!*******************************************************************************
! converts from lat/long to National Grid numbers 
! East,North are given in metres
! Lat,Long are required, as:
! EITHER   in degrees as a decimal, with minutes and seconds both set to zero
! OR       in degrees, minutes, seconds

subroutine LatLongToNatGrid (East,North,LatDeg,LatMin,LatSec,LonDeg,LonMin,LonSec)

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

real (kind=big), intent (out) :: East,North
real (kind=big), intent (in)  :: LatDeg,LatMin,LatSec,LonDeg,LonMin,LonSec

real (kind=big), parameter :: Ka = 6377563.396
real (kind=big), parameter :: Kb = 6356256.910
real (kind=big), parameter :: Kpi = 3.14159265358979
real (kind=big), parameter :: KN0 = -100000.0
real (kind=big), parameter :: KE0 =  400000.0
real (kind=big), parameter :: KF0 = 0.9996012717
real (kind=big), parameter :: KThi0 = 49.0
real (kind=big), parameter :: KLam0 = -2.0

real (kind=big) :: Lat, Lon
real (kind=big) :: Ve2, Vn,Vv, Vro, Vnu2, VM, V1, V2, V3, V3A, V4, V5, V6, V6A, V6B, V6C, Br, Br2,VThi0,VLam0

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

Lat = LatDeg + (LatMin/real(60,big)) + (LatSec/real(3600,big))
Lon = LonDeg + (LonMin/real(60,big)) + (LonSec/real(3600,big))

Lat = (Lat * Kpi) / real(180,big)
Lon = (Lon * Kpi) / real(180,big)

VThi0 = (KThi0 / real(180,big)) * Kpi
VLam0 = (KLam0 / real(180,big)) * Kpi
Ve2   = ((Ka*Ka)-(Kb*Kb))/(Ka*Ka)
Vn    = (Ka-Kb) / (Ka+Kb)

Vv   = Ka * KF0 * ((real(1,big) - (Ve2 * sin(Lat) * sin(Lat))) ** (0-0.5))
Vro  = Ka * KF0 * (real(1,big) - Ve2) * ((real(1,big) - (Ve2 * sin(Lat) * sin(Lat))) ** (0-1.5))
Vnu2 = (Vv/Vro) - real(1,big)

VM = CalcM (Ka,Kb,KF0,VThi0,Lat)

V1  = VM + KN0
V2  = (Vv/real(2,big))*sin(Lat)*cos(Lat)
V3  = (Vv/real(24,big))*sin(Lat)*(cos(Lat)**3)*(real(5,big)-(tan(Lat)*tan(Lat))+(real(9,big)*Vnu2))
V3A = (Vv/real(720,big))*sin(Lat)*(cos(Lat)**5)*(real(61,big)-(real(58,big)*(tan(Lat)**2))+(tan(Lat)**4))
V4  = (Vv*cos(Lat))
V5  = (Vv/real(6,big))*(cos(Lat)**3)*((Vv/Vro)-(tan(Lat)**2))
V6A = (Vv/real(120,big))*(cos(Lat)**5)
V6B = real(5,big)-(real(18,big)*(tan(Lat)**2))+(tan(Lat)**4)+(real(14,big)*Vnu2)
V6C = real(58,big)*(tan(Lat)**2)*Vnu2
V6  = V6A * (V6B - V6C)

Br  = Lon - VLam0
Br2 = Br * Br

North = V1 + (V2*Br2) + (V3*Br2*Br2) + (V3A*Br2*Br2*Br2)
East  = KE0 + (V4*Br) + (V5*Br*Br2) + (V6*Br*Br2*Br2)

end subroutine LatLongToNatGrid

!*******************************************************************************
! calcs 'M'

function CalcM (Va,Vb,VF0,VThi0,VThi)

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

real (kind=big) :: CalcM,Vn
real (kind=big), intent (in) :: Va,Vb,VF0,VThi0,VThi

Vn = (Va-Vb) / (Va+Vb)

CalcM =         (real(1,big)+Vn+(real(1.25,big)*Vn*Vn)+(real(1.25,big)*Vn*Vn*Vn)) * (VThi-VThi0)

CalcM = CalcM - (((real(3,big)*Vn)+(real(3,big)*Vn*Vn)+(real((2.625),big)*Vn*Vn*Vn)) * sin(VThi-VThi0) * cos(VThi+VThi0))

CalcM = CalcM + (((real(1.875,big)*Vn*Vn)+(real(1.875,big)*Vn*Vn*Vn)) * sin(real(2,big)*(VThi-VThi0)) * cos(real(2,big)*(VThi+VThi0)))

CalcM = CalcM - (((real(35,big)/real(24,big))*Vn*Vn*Vn) * sin(real(3,big)*(VThi-VThi0)) * cos(real(3,big)*(VThi+VThi0)))

CalcM = CalcM * Vb * VF0

end function CalcM

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

end module GridOps