I Overview ---------- Computation of slant path delay through the neutral atmosphere for reduction of SLR, GPS, VLBI and DORIS data is performed by three procedures. The first procedure, data acquisition is performed by the designated server. The server downloads the numerical weather model GEOS-FPIT (Rienecker et al. 2008) developed and maintained by the Goddard Modeling and Assimilation Office. The model has resolution 0.625 deg x 0.5 deg x 72 layers x 3 hours, runs from 2000.01.01 through presents, updated 4 times a day and has latency 10-15 hours. The output of the numerical weather model among other parameters contains thickness of atmospheric layers, air temperature, and specific humidity at each grid point. These parameters are used for computing refractivities at three wavelength ranges: 532nm, 1064 nm, and 1mm-30m (radio). The refractivity accounts for both dry air and water vapor. The data acquisition procedure stores for each epoch 3D refractivity field at a regular grid with the height range from -1000, to 80000 km above the reference ellipsoid. The data acquisition procedure runs every hour and check whether the new output from the GEOS-FPIT model become available. The second procedure computes slant path delay for a designated set of stations. The path delay are computed on a regular, but non-uniform grid over azimuth and elevation. For each station, each epoch, each direction the trajectory of the wavefront is computed by solving a system of non-liner differential equations of the 4th order that are the solution of the variational problem of wave propagation in accordance with the Fermat principle. Then the path delay is computed by integrating the refractivity along the trajectory from the receiver (i.e. observing station) to the top of the atmosphere defined at the height of 80 km. All slant path delay for all directions, all stations for a given epochs are stored in an output file. Path delays for designated "continuous" stations are computed immediately after completion of the data acquisition procedure. The third procedure is incorporated into a space geodesy data reduction software, such as GEODYN and Calc/Solve. The procedure determines the time range of observations and the list of stations. It downloads files with slant path delay on an azimuth-elevation grid. For each station it expands the path delay over the tensor product of the 3D B-spline basis that runs over azimuth, elevation, and time. Then, using these expansion coefficients, slant path delay and its partial derivatives with respect to the path delay in zenith direction are computed to a given moment of time, given elevation and azimuth of the geodetic observation. The essential pre-requisite in this model is that the position of the receiver is known with accuracy 1 meter and the emitter is assumed to be at distance greater than 1000 km. II. Routines: The package SPD_CLIENT provides two routines for the the third procedure that are supposed to be called from a space geodesy data reduction software. The first routine, SPD_LOAD_BSPD, loads the data from a directory with expansion coefficients in bspd format. The second routine, SPD_INTRP_DELAY, compute path delay by interpolation. a) SPD_LOAD_BSPD: Routine SPD_LOAD_BSPD reads a list of directories with slant path delays in binary format, finds files for the specified list of stations ad the specified date range [MJD_BEG/TAI_BEG, MJD_END/TAI_END]. It reads the path delays and for each station computes te coefficients of the 3D B-spline expansion of path delay over mapping value, azimuth and time. Mapping values is a function of elevation angle. It is the ratio of path delay at a given elevation angle to the path delay at zenith direction for the ISA 1976 Standard Atmosphere. The coefficients are written in fields of array of objects SPD. Slant path delay for a given station may be present in several directories. If so, the definition in the first directories will take the precedence. That mean, if the data for a given station present in i-th directory, the i+1, i+2, etc will not be sought. ________________________ Input parameters: _________________________ L_DIR ( INTEGER*4 ) -- The number of directories with slant path delays in BSPD format. C_DIR ( CHARACTER ) -- List of directories with slant path delay in BSPD format. Dimension: L_DIR. SPD_BIAS_FILE ( CHARACTER ) -- File with biases of slant path delay. NONE means no biases to slant path delay will be applied. APD_NAME ( CHARACTER ) -- Model of partial derivatives with respect to atmosphere path delay in zenith direction. Supported models: NMFW -- Niell (1996) mapping function (wet). NMFH -- Niell (1996) mapping function (hydrostatic) on the mean epoch. TOTAL_SCALE -- mapping function for the case when residual atmosphere is considered proportional to the total atmosphere. The partial derivative is defined as a ratio of the total slant path delay to the total path delay in the zenith direction. WATER_SCALE -- mapping function for the case when residual atmosphere is considered proportional to the water vapor contribution to the atmosphere. The partial derivative is defined as a ratio of the water vapor contribution of slant path delay to the water vapor contribution to path delay in the zenith direction. GAUSSIAN_LAYER -- mapping function for the case when the dependence of concentration of the residual atmosphere is described with the Gaussian model with the specified height and the specified full width half maximum (FWHM). LAYER_HEIGHT ( REAL*8 ) -- If the partial derivative model is GAUSSIAN_LAYER, then APD_PAR1 is the layer height in meters. Otherwise, it is zero. LAYER_FWHM ( REAL*8 ) -- If the partial derivative model is GAUSSIAN_LAYER, then APD_PAR1 is the layer FWHM in meters. Otherwise, it is zero. L_STA ( INTEGER*4 ) -- The number of stations for which slant path delay will be interpolated. C_STA ( CHARACTER ) -- A list of station names for which slant path delay will be interpolated. The name should be in the list of stations for which the path delay have been computed. Dimension: L_STA. MJD_BEG ( INTEGER*4 ) -- Integer Modified Julian Date of the beginning of the interval of interpolation. Units: days. TAI_BEG ( REAL*8 ) -- TAI time since the midnight of the beginning of the interval of interpolation. Units: sec. MJD_END ( INTEGER*4 ) -- Integer Modified Julian Date of the end of the interval of interpolation. Units: days. TAI_END ( REAL*8 ) -- TAI time since the midnight of the end of the interval of interpolation. Units: sec. _________________________ Output parameters: _______________________ SPD ( SPD_DEL__TYPE ) -- Array of objects with Slant Path Delays. Routine SPD_LOAD_BSPD will populate internal fields of SPD array. Dimension: L_STA. ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. b) SPD_LOAD_BSPD: Routine SPD_INTRP_DELAY computes slant path delay, slant path delay rate, partial derivative of slant path delay on the delay in zenith direction, and partial derivative of slant path delay on the delay in zenith direction. These computations are based on the coefficients of expansion of path delay over 3D B-spline basis. It is assumed routine SPD_LOAD_BSPD was been called before, it has computed the coefficients of expansion slant path delay over B-spline basis and written them into array of SPD objects for a number of stations, including the station STA_NAM in this routine. ________________________ Input parameters: _________________________ STA_NAM ( CHARACTER ) -- Station name. The station should be from the list supplied to SPD_LOAD_BSPD. L_STA ( INTEGER*4 ) -- The number of stations for which the B-spline expansion coefficients has been computed by routine SPD_LOAD_BSPD. SPD ( SPD_DEL__TYPE ) -- Array of objects with Slant Path Delays. The internal fields of SPD array contains the coefficient of B-spline expansion of path delay. Dimension: L_STA. ELEV ( REAL*8 ) -- Elevation angle. Units: radian. AZIM ( REAL*8 ) -- Azimuth counted from North to East. Units: radian. ELEV_RATE ( REAL*8 ) -- Rate of changes of elevation angle. Units: rad/s. AZIM_RATE ( REAL*8 ) -- Rate of changes of azimuth angle. Units: rad/s. MJD_OBS ( INTEGER*4 ) -- Integer Modified Julian Date of the observation at midnight. Units: days. TAI_BEG ( INTEGER*4 ) -- TAI time since the midnight of the observation. Units: sec. _________________________ Output parameters: _______________________ SPD_DELAY ( REAL*8 ) -- Slant path delay. Units: s. SPD_RATE ( REAL*8 ) -- Time derivative of slant path delay. Units: dimensionless SPD_DELAY_DER_ZEN ( REAL*8 ) -- Partial derivative of slant path delay at a given AZ,EL over the path delay in zenith direction using only wet contribution to path delay (mapping function). Units: dimensionless. SPD_RATE_DER_ZEN ( REAL*8 ) -- Partial derivative of slant path rate at a given AZ,EL over the path delay in zenith direction using only wet contribution to path delay. Units: 1/sec ________________________ Modified parameters: ______________________ IUER ( INTEGER*4, OPT ) -- Universal error handler. Input: switch IUER=0 -- no error messages will be generated even in the case of error. IUER=-1 -- in the case of error the message will be put on stdout. Output: 0 in the case of successful completion and non-zero in the case of error. 3) Auxilary program spd_3d_toser The second procedure computes path delay in a file in ascii spd_3d format , one file per epoch. That file contains path delay for all the stations of a given epoch. SPD_LOAD_BSPD expects to find slant path in a file in binary format bspd. The data in bspd formats holds all path delays for a given station. The server that computes path delays in ascii format also also invokes program spd_3d_toser that transforms the results in binary bspd format. That program is also provided by the SPD_CLIENT package. Usage: spd_3d_toser {dir_in} {pref_out} {create/update} [ivrb] 1st argument dir_in is the directory with files with path delays in spd_3d format. The second argument is the prefix that is prepended to the output files, usually full path directory name and some word. The third argument is either create or update. If create is selected the new output files will be created. If update is selected, the output files will be read and new records will be appended to the end. The fourth optional argument is verbosity level. irvb=0 means silent mode, ivrb=1 means normal verbosity (default), and ivrb=2 shows the progress of the operation. NB: No gaps in the time series of slant path delay is allowed. If a gap will be found, an error message will be issued. III. Example: ------------- spd_show -- example of computing path delay at a given epoch. IV: References: --------------- Petrov, L., ICESAT-2 Algorithm Technical Base Document for Atmospheric delay correction to laser altimetry ranges, 2015, submitted Petrov, L., Modeling of path delay in the neutral atmosphere: a paradigm shift, to appear in the Proceedings of the 12th European VLBI Network Symposium and Users Meeting, 7-10 October 2014 Cagliari, Italy http://arxiv.org/abs/1502.06678 Rienecker, M.M., Suarez, M.J., Todling, R., Bacmeister, J., Takacs, L., Liu, H.C., Sienkiewicz, W.M., Koster, R.D., Gelaro, R., Stajner I., and Nielsen E., ``The GEOS Data Assimilation System --- Documentation of Versions 5.0.1, 5.1.0, and 5.2.0.'', NASA/TM--2008--104606, 2008. http://gmao.gsfc.nasa.gov/pubs/docs/tm27.pdf