GeodeticShiftToWGS84 Subroutine

private subroutine GeodeticShiftToWGS84(input, latIn, lonIn, WGS84lat, WGS84lon)

shifts geodetic coordinates relative to a given source datum to geodetic coordinates relative to WGS84.

Arguments

Type IntentOptional Attributes Name
type(Coordinate), intent(in) :: input
real(kind=float), intent(in) :: latIn
real(kind=float), intent(in) :: lonIn
real(kind=float), intent(out) :: WGS84lat
real(kind=float), intent(out) :: WGS84lon

Variables

Type Visibility Attributes Name Initial
real(kind=float), public, parameter :: MOLODENSKY_MAX = 89.75*degToRad

Polar limit

real(kind=float), public :: WGS84_a

Semi-major axis of WGS84 ellipsoid in meters

real(kind=float), public :: WGS84_f

Flattening of WGS84 ellisoid

real(kind=float), public :: WGS84hgt
real(kind=float), public :: a

Semi-major axis of ellipsoid in meters

real(kind=float), public :: da

Difference in semi-major axes

real(kind=float), public :: df

Difference in flattening

real(kind=float), public :: dx
real(kind=float), public :: dy
real(kind=float), public :: dz
real(kind=float), public :: f

Flattening of ellipsoid

real(kind=float), public :: hgtIn

Source Code

SUBROUTINE GeodeticShiftToWGS84 &
!
(input, latIn, lonIn, WGS84lat, WGS84lon)

IMPLICIT NONE

! Arguments with intent(in):
TYPE (Coordinate), INTENT(IN) :: input
REAL (KIND = float), INTENT(IN) :: latIn
REAL (KIND = float), INTENT(IN) :: lonIn

! Arguments with intent(out):
REAL (KIND = float), INTENT(OUT) :: WGS84lat, WGS84lon

! Local variables:
REAL (KIND = float) :: hgtIn
REAL (KIND = float) :: WGS84_a !! Semi-major axis of WGS84 ellipsoid in meters 
REAL (KIND = float) :: WGS84_f !! Flattening of WGS84 ellisoid    
REAL (KIND = float) :: a !!Semi-major axis of ellipsoid in meters
REAL (KIND = float) :: f !!Flattening of ellipsoid
REAL (KIND = float) :: da !!Difference in semi-major axes
REAL (KIND = float) :: df !!Difference in flattening
REAL (KIND = float) :: dx
REAL (KIND = float) :: dy
REAL (KIND = float) :: dz 
REAL (KIND = float) :: WGS84hgt


! Local parameters:
REAL (KIND = float), PARAMETER :: MOLODENSKY_MAX  = 89.75 * degToRad !! Polar limit          

!------------end of declaration------------------------------------------------

!initialize
hgtIn = 0.
WGS84_a = 6378137.0 
WGS84_f = 1. / 298.257223563
a = input % system % ellipsoid % a
f = input % system % ellipsoid % f

!If input datum is WGS84, just copy to output
IF ( input % system % datum % code == WGS84 ) THEN
  WGS84lat = latIn
  WGS84lon = lonIn
  RETURN
END IF

!If latitude is within limits, apply Molodensky 
IF ( latIn <= MOLODENSKY_MAX .AND. &
     latIn >= - MOLODENSKY_MAX  ) THEN
     
   da = WGS84_a - a
   df = WGS84_f - f
   dx = input % system % datum % dx
   dy = input % system % datum % dy
   dz = input % system % datum % dz
   CALL MolodenskyShift(WGS84_a, da, WGS84_f, df, dx, dy, dz, latIn, lonIn, &
                         hgtIn, WGS84lat, WGS84lon, WGS84hgt)
     
ELSE !apply 3 steps method

  WRITE (*,*) 'Latitude is outside Molodensky limits'
  WRITE (*,*) '3-step datum transformation not yet implemented'
  STOP

END IF

END SUBROUTINE GeodeticShiftToWGS84