subroutine saasthyd (p,dlat,hell,zhd)

c% This subroutine determines the zenith hydrostatic delay based on the
c% equation by Saastamoinen (1972) as refined by Davis et al. (1985)
c%
c% c Reference:
c% Saastamoinen, J., Atmospheric correction for the troposphere and 
c% stratosphere in radio ranging of satellites. The use of artificial 
c% satellites for geodesy, Geophys. Monogr. Ser. 15, Amer. Geophys. Union, 
c% pp. 274-251, 1972.
c% Davis, J.L, T.A. Herring, I.I. Shapiro, A.E.E. Rogers, and G. Elgered, 
c% Geodesy by Radio Interferometry: Effects of Atmospheric Modeling Errors 
c% on Estimates of Baseline Length, Radio Science, Vol. 20, No. 6, 
c% pp. 1593-1607, 1985.
c%
c% input parameters:
c%
c% p:     pressure in hPa
c% dlat:  ellipsoidal latitude in radians 
c% dlon:  longitude in radians 
c% hell:  ellipsoidal height in m 
c% 
c% output parameters:
c%
c% zhd:  zenith hydrostatic delay in meter 
c%
c% Example 1 :
c%
c% p = 1000;
c% dlat = 48d0*pi/180.d0
c% hell = 200.d0
c%
c% output:
c% zhd = 2.2695 m
c%
c% Johannes Boehm, 8 May 2013
c% Johannes Boehm, 24 December 2014, converted to Fortran
c% ---

      implicit double precision (a-h,o-z)

      !% calculate denominator f
      f = 1.d0-0.00266d0*dcos(2.d0*dlat) - 0.00000028d0*hell

      !% calculate the zenith hydrostatic delay
      zhd = 0.0022768d0*p/f

      end subroutine