<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
!/**----------------------------------------------------------------------
! @sub msl2geo1
!
! Convert a list of geometric mean sea level altitudes to geopotential height
! Use a more accurate conversion found at:
! http://mtp.jpl.nasa.gov/notes/altitude/altitude.html
!
! @parameter
! @ input: z -- list of geometric altitudes (in km)
! @ lat -- latitude of profile in degrees
! @ lon -- longitude of profile in degrees.
! @ output: h -- list of geopotential altitudes, in km
! -----------------------------------------------------------------------*/
<A NAME='DA_MSL2GEO1'><A href='../../html_code/tools/da_msl2geo1.inc.html#DA_MSL2GEO1' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_msl2geo1 (z, lat, lon, h) 1
implicit none
real*8 h, lat, lon, z ! Note lon is currently not used
real*8 pi, latr, zm
real*8 semi_major_axis, semi_minor_axis, grav_polar, grav_equator
real*8 earth_omega, grav_constant, flattening, somigliana
real*8 grav_ratio, sin2, termg, termr, grav, eccentricity
! Parameters below from WGS-84 model software inside GPS receivers.
parameter(semi_major_axis = 6378.1370d3) ! (m)
parameter(semi_minor_axis = 6356.7523142d3) ! (m)
parameter(grav_polar = 9.8321849378) ! (m/s2)
parameter(grav_equator = 9.7803253359) ! (m/s2)
parameter(earth_omega = 7.292115d-5) ! (rad/s)
parameter(grav_constant = 3.986004418d14) ! (m3/s2)
parameter(grav = 9.80665d0) ! (m/s2) WMO std g at 45 deg lat
parameter(eccentricity = 0.081819d0) ! unitless
! Derived geophysical constants
parameter(flattening = (semi_major_axis-semi_minor_axis) / &
semi_major_axis)
parameter(somigliana = &
(semi_minor_axis/semi_major_axis)*(grav_polar/grav_equator)-1.d0)
parameter(grav_ratio = (earth_omega*earth_omega * &
semi_major_axis*semi_major_axis * semi_minor_axis)/grav_constant)
pi = 3.14159265358979d0
latr = lat * (pi/180.d0) ! in radians
zm = z*1000d0 ! in meters
if (z.eq.-999.d0) then
h = -999.d0
else
sin2 = sin(latr) * sin(latr)
termg = grav_equator * &
( (1.d0+somigliana*sin2) / &
sqrt(1.d0-eccentricity*eccentricity*sin2) )
termr = semi_major_axis / &
(1.d0 + flattening + grav_ratio - 2.d0*flattening*sin2)
! geopotential height
h=((termg/grav)*((termr*zm)/(termr+zm)))*0.001 ! in km
endif
end subroutine da_msl2geo1