earth_radius.f90

 SUBROUTINE earth_radius(IGDTMPL, IGDTLEN, RADIUS, ECCEN_SQUARED)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM:  EARTH_RADIUS   DETERMINE EARTH RADIUS AND SHAPE
!   PRGMMR: GAYNO    ORG: W/NMC23     DATE: 2015-07-14
!
! ABSTRACT: DETERMINE THE RADIUS AND SHAPE OF THE EARTH FROM
!   THE GRIB 2 GRID DEFINITION TEMPLATE ARRAY - SECTION 3
!
! PROGRAM HISTORY LOG:
! 2015-07-14  GAYNO
!
! USAGE:   CALL EARTH_RADIUS(IGDTMPL, IGDTLEN, &
!                            RADIUS, ECCEN_SQUARED)
!
!   INPUT ARGUMENT LIST:
!     IGDTMPL  - INTEGER (IGDTLEN) GRID DEFINITION TEMPLATE ARRAY.
!                CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT
!                OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!                FOR ALL MAP PROJECTIONS RECOGNIZED BY IPLIB,
!                THE ENTRIES USE BY THIS ROUTINE ARE:
!                 (1) - SHAPE OF EARTH, SECTION 3, OCTET 15 
!                 (2) - SCALE FACTOR OF SPHERICAL EARTH RADIUS,
!                       OCTET 16
!                 (3) - SCALED VALUE OF RADIUS OF SPHERICAL EARTH,
!                       OCTETS 17-20
!                 (4) - SCALE FACTOR OF MAJOR AXIS OF ELLIPTICAL EARTH,
!                       OCTET 21
!                 (5) - SCALED VALUE OF MAJOR AXIS OF ELLIPTICAL EARTH,
!                       OCTETS 22-25
!                 (6) - SCALE FACTOR OF MINOR AXIS OF ELLIPTICAL EARTH,
!                       OCTET 26
!                 (7) - SCALED VALUE OF MINOR AXIS OF ELLIPTICAL EARTH,
!                       OCTETS 27-30
!     IGDTLEN  - INTEGER NUMBER OF ELEMENTS OF THE GRID DEFINITION
!                TEMPLATE ARRAY.  CORRESPONDS TO THE GFLD%IGDTLEN
!                COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD 
!                DATA STRUCTURE.
!
!   OUTPUT ARGUMENT LIST:
!     RADIUS        - REAL EARTH RADIUS IN METERS.
!                     FOR ELLIPITICAL EARTHS, THIS IS THE
!                     SEMI MAJOR AXIS.  SEE "MAP PROJECTSIONS -
!                     A WORKING MANUAL" BY SNYDER (1987)
!                     FOR DETAILS.
!     ECCEN_SQUARED - REAL EARTH ECCENTRICITY SQUARED
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
 IMPLICIT NONE
!
 INTEGER,                INTENT(IN   ) :: IGDTLEN
 INTEGER,                INTENT(IN   ) :: IGDTMPL(IGDTLEN)
!
 REAL,                   INTENT(  OUT) :: ECCEN_SQUARED
 REAL,                   INTENT(  OUT) :: RADIUS
!
 REAL                                  :: FLAT
 REAL                                  :: MAJOR_AXIS, MINOR_AXIS
!
 SELECT CASE (igdtmpl(1))
   CASE (0)
     radius        = 6367470.0
     eccen_squared = 0.0
   CASE (1)  ! USER SPECIFIED SPHERICAL
     radius        = float(igdtmpl(3))/float(10**igdtmpl(2))
     eccen_squared = 0.0
   CASE (2)  ! IAU 1965
     radius        = 6378160.0      ! SEMI MAJOR AXIS
     flat          = 1.0/297.0      ! FLATTENING
     eccen_squared = (2.0*flat) - (flat**2)
   CASE (3)  ! USER SPECIFIED ELLIPTICAL (KM)
     major_axis    = float(igdtmpl(5))/float(10**igdtmpl(4))
     major_axis    = major_axis * 1000.0
     minor_axis    = float(igdtmpl(7))/float(10**igdtmpl(6))
     minor_axis    = minor_axis * 1000.0
     eccen_squared = 1.0 - (minor_axis**2 / major_axis**2)
     radius        = major_axis
   CASE (4)  ! IAG-GRS80 MODEL
     radius        = 6378137.0      ! SEMI MAJOR AXIS
     flat          = 1.0/298.2572   ! FLATTENING
     eccen_squared = (2.0*flat) - (flat**2)
   CASE (5)  ! WGS84 DATUM
     radius        = 6378137.0      ! SEMI MAJOR AXIS
     eccen_squared = 0.00669437999013
   CASE (6)
     radius        = 6371229.0
     eccen_squared = 0.0
   CASE (7)  ! USER SPECIFIED ELLIPTICAL (M)
     major_axis    = float(igdtmpl(5))/float(10**igdtmpl(4))
     minor_axis    = float(igdtmpl(7))/float(10**igdtmpl(6))
     eccen_squared = 1.0 - (minor_axis**2 / major_axis**2)
     radius        = major_axis
   CASE (8)
     radius        = 6371200.0
     eccen_squared = 0.0
   CASE DEFAULT
     radius        = -9999.
     eccen_squared = -9999.
 END SELECT
!
 RETURN
!
 END SUBROUTINE earth_radius