var/da/da_tools/da_geo2msl1.inc

   1 !/**----------------------------------------------------------------------
   2 ! @sub       geo2msl1
   3 !
   4 ! Convert a geopotential height to mean sea level geometric height.
   5 ! Use a more accurate conversion found at:
   6 ! http://mtp.jpl.nasa.gov/notes/altitude/altitude.html
   7 !
   8 ! @parameter
   9 ! @ input:    h   -- geopotential height (km)
  10 ! @           lat -- latitude  of profile in degrees
  11 ! @           lon -- longitude of profile in degrees.
  12 ! @ output:   z   -- geometric height (km)
  13 ! -----------------------------------------------------------------------*/
  14 subroutine da_geo2msl1 (h, lat, lon, z)
  15       implicit none
  16       real*8 h, lat, lon, z     ! Note lon is currently not used
  17       real*8 pi, latr, hm
  18       real*8 semi_major_axis, semi_minor_axis, grav_polar, grav_equator
  19       real*8 earth_omega, grav_constant, flattening, somigliana
  20       real*8 grav_ratio, sin2, termg, termr, grav, eccentricity
  21
  22 !     Parameters below from WGS-84 model software inside GPS receivers.
  23       parameter(semi_major_axis = 6378.1370d3)    ! (m)
  24       parameter(semi_minor_axis = 6356.7523142d3) ! (m)
  25       parameter(grav_polar = 9.8321849378)        ! (m/s2)
  26       parameter(grav_equator = 9.7803253359)      ! (m/s2)
  27       parameter(earth_omega = 7.292115d-5)        ! (rad/s)
  28       parameter(grav_constant = 3.986004418d14)   ! (m3/s2)
  29       parameter(grav = 9.80665d0)                 ! (m/s2) WMO std g at 45 deg lat
  30       parameter(eccentricity = 0.081819d0)        ! unitless
  31
  32 !     Derived geophysical constants
  33       parameter(flattening = (semi_major_axis-semi_minor_axis) /  &
  34                               semi_major_axis)
  35
  36       parameter(somigliana = &
  37        (semi_minor_axis/semi_major_axis)*(grav_polar/grav_equator)-1.d0)
  38
  39       parameter(grav_ratio = (earth_omega*earth_omega * &
  40        semi_major_axis*semi_major_axis * semi_minor_axis)/grav_constant)
  41
  42       pi   = 3.14159265358979d0
  43       latr = lat * (pi/180.d0)        ! in radians
  44       hm   = h * 1000d0               ! in meters
  45
  46       if (h.eq.-999.d0) then
  47          z = -999.d0
  48       else
  49          sin2  = sin(latr) * sin(latr)
  50          termg = grav_equator *  &
  51           ( (1.d0+somigliana*sin2) /  &
  52             sqrt(1.d0-eccentricity*eccentricity*sin2) )
  53          termr = semi_major_axis   /  &
  54           (1.d0 + flattening + grav_ratio - 2.d0*flattening*sin2)
  55
  56 !     geometric height
  57
  58          z = ( (termr*hm) / ( (termg/grav) * termr - hm ) ) * 0.001 ! in km
  59
  60       endif
  61
  62 end subroutine da_geo2msl1