ConvertGeodeticToUTM Subroutine

SUBROUTINE ConvertGeodeticToUTM &
!
(lon, lat, k, centM, lat0, a, e, eb, falseN, falseE, x, y, zone, &
 hemisphere, override)

USE Units, ONLY: &
!Imported parametes:
pi

USE StringManipulation, ONLY: &
!Imported routines:
ToString

IMPLICIT NONE

! Arguments with intent(in):
REAL (KIND = float), INTENT (IN) :: k !!scale factor
REAL (KIND = float), INTENT (IN) :: lat0 !!latitude of origin [radians]
REAL (KIND = float), INTENT (IN) :: a !! semimajor axis [m]
REAL (KIND = float), INTENT (IN) :: e !! eccentricity
REAL (KIND = float), INTENT (IN) :: eb !! second eccentricity
REAL (KIND = float), INTENT (IN) :: falseE !!false easting
REAL (KIND = float), INTENT (IN) :: override !!override option

!Arguments with intent(inout):
REAL (KIND = float), INTENT (INOUT) :: lon !!geodetic longitude [radians]
REAL (KIND = float), INTENT (INOUT) :: lat !!geodetic latitude [radians]
REAL (KIND = float), INTENT (INOUT) :: centM !!central meridian [radians]
REAL (KIND = float), INTENT (INOUT) :: falseN !!false northing
REAL (KIND = float), INTENT (INOUT) :: zone !!zone recalculated if override = 0.
REAL (KIND = float), INTENT (INOUT) :: hemisphere !!hemisphere recalculated if override = 0.

!Arguments with intent(out):
REAL (KIND = float), INTENT (OUT) :: x !!easting coordinate [m]
REAL (KIND = float), INTENT (OUT) :: y !!northing coordinate [m]

!Local variables:
INTEGER (KIND = long) :: lat_degrees
INTEGER (KIND = long) :: lon_degrees

!Local parameters:
REAL (KIND = float), PARAMETER :: MIN_LAT = -80.5 * degToRad ! -80.5 degrees in radians 
REAL (KIND = float), PARAMETER :: MAX_LAT = 84.5 * degToRad ! 84.5 degrees in radians 
REAL (KIND = float), PARAMETER :: MIN_EASTING = 100000.
REAL (KIND = float), PARAMETER :: MAX_EASTING = 900000.
REAL (KIND = float), PARAMETER :: MIN_NORTHING = 0.
REAL (KIND = float), PARAMETER :: MAX_NORTHING = 10000000.
!------------end of declaration------------------------------------------------
!check out of range
IF ( lat < MIN_LAT .OR. lat > MAX_LAT ) THEN
  CALL Catch ('error', 'GeoLib',   &
			 'Converting Geodetic to UTM: &
			 latitude out of range ' ,  &
			 code = consistencyError, &
			 argument = ToString(lat*radToDeg)//' deg' )
END IF

IF ( lon < -pi .OR. lon > 2*pi ) THEN
  CALL Catch ('error', 'GeoLib',   &
			 'Converting Geodetic to UTM: &
			 longitude out of range ' ,  &
			 code = consistencyError, &
			 argument = ToString(lon*radToDeg)//' deg' )
END IF

IF ( lat > -1.0e-9 .AND. lat < 0. ) THEN
  lat = 0.0
END IF

IF ( lon < 0. ) THEN
  lon = lon + 2.*pi + 1.0e-10
END IF

!calculate latitude and longitude in degrees
lat_degrees = INT (lat * radToDeg)
lon_degrees = INT (lon * radToDeg)

!If override is off, recalculate zone and centM
IF ( override == 0. ) THEN
  IF (lon < pi) THEN
    zone = INT ( (31. + ((lon * radToDeg) / 6.0)) )
  ELSE
    zone = INT ( ( (lon * radToDeg) / 6.0 - 29.) )
  END IF
  
  IF (zone > 60.) THEN
      zone = 1.
  END IF
  
  !UTM special cases 
  IF ( (lat_degrees > 55) .AND. (lat_degrees < 64) .AND. &
       (lon_degrees > -1) .AND. (lon_degrees < 3) ) THEN
      zone = 31.
  END IF
  IF ( (lat_degrees > 55) .AND. (lat_degrees < 64) .AND. &
       (lon_degrees > 2) .AND. (lon_degrees < 12) ) THEN
      zone = 32.
  END IF
  IF ( (lat_degrees > 71) .AND. (lon_degrees > -1) .AND. &
       (lon_degrees < 9) ) THEN
      zone = 31.
  END IF
  IF ( (lat_degrees > 71) .AND. (lon_degrees > 8) .AND. &
       (lon_degrees < 21) ) THEN
      zone = 33.
  END IF
  IF ( (lat_degrees > 71) .AND. (lon_degrees > 20) .AND. &
       (lon_degrees < 33) ) THEN
      zone = 35.
  END IF
  IF ( (lat_degrees > 71) .AND. (lon_degrees > 32) .AND. &
       (lon_degrees < 42) ) THEN
      zone = 37.
  END IF
          
END IF

IF (zone >= 31.) THEN
   centM = (6. * zone - 183.) * degToRad
ELSE
   centM = (6. * zone + 177.) * degToRad
END IF

IF (lat < 0.) THEN
  falseN = 10000000.
  hemisphere = SOUTH
ELSE
  falseN = 0.
  hemisphere = NORTH
END IF


CALL ConvertGeodeticToTransverseMercator (lon, lat, k, centM, lat0, a, e, eb, &
                                          falseN, falseE, x, y)
                                          
!Check out of range if override is off
IF ( override == 0. ) THEN
    IF ( x < MIN_EASTING .OR. x > MAX_EASTING ) THEN
      CALL Catch ('error', 'GeoLib',   &
			      'Converting Geodetic to UTM: &
			       easting out of range' ,  &
			       code = consistencyError, argument = ToString(x) )
    END IF
    IF ( y < MIN_NORTHING .OR. y > MAX_NORTHING ) THEN
      CALL Catch ('error', 'GeoLib',   &
			      'Converting Geodetic to UTM: &
			       northing out of range' ,  &
			       code = consistencyError, argument = ToString(y) )
    END IF
END IF
      
END SUBROUTINE ConvertGeodeticToUTM