Ray tracing at CIRES/NOAA
(as of 13 July 2018)

To the left is an example of the effect of wind and temperature structure on the propagation of acoustic waves in the atmosphere. We have developed computer programs to calculate the propagation of acoustic waves (with or without winds and attenuation loss) and acoustic-gravity waves (without attenuation loss and with or without wind) in the atmosphere or ocean.

In addition, we have programs to plot contours and/or profiles of sound speed, temperature, etc. for the corresponding atmospheric or oceanic models. All of these programs are available for downloading.

Other examples of some raypaths we have calculated are given here and here. More information can be found in the HARPA manual, the HARPO manual, and our ray tracing publications. A tutorial on ray tracing can be found in "Three-dimensional ray tracing in the atmosphere." , which is reproduced here. The errata are at http://cires.colorado.edu/~mjones/pubs/errata2.htm .

A tutorial on dispersion relations is here. Further discussion of the propagation of gravity waves can be found in our publications on internal gravity-wave propagation. A dispersion-relation diagram is sometimes referred to as a propagation-surface diagram.

An explanation of the coordinate systems used to plot raypaths and make contour plots is given here. An explanation of the machine-readable output is given here, both for the transmitter raysets and the receiver raysets .

A list of the subroutine and function names is given here.

CONTENTS

  1. Instructions for downloading the ray tracing program .
  2. Directories for the ray tracing program and supporting programs .
  3. Earlier versions of some of our programs .
  4. Notice about numerical integration errors and computer word length .
  5. How to find out what the integration error is .
  6. The gravity-wave dispersion relation with wind .
  7. Dispersion relation for magneto-acoustic-gravity waves .
  8. Algorithms we use to find intersections of the ray with the terrain and with the receiver surface .
  9. Guidelines we have tried to follow in developing our ray tracing programs .
  10. Updates to the ray tracing programs .
  11. New models in subdirectory /raytracing/models/ .
  12. Profile and contour plotting .
  13. Improvements still to be made to the ray tracing program .
  14. Improvements to the documentation of the ray tracing program .
  15. How to add new models to the ray tracing program .

The raytracing directory has all of the ray tracing programs that we are distributing. Each program is in a separate subdirectory. Most of those subdirectories have detailed descriptions. Here, we give brief descriptions of the programs we have written.

If you want to run these programs, here is what you need to do.

  1. Download all of the following directories with all of their files.
  2. Compile (but do not link) all of the .for files in the directories pcharpo, models, conplt, and g.
  3. Compile (and link) eigen.for (using the batch file elnk.bat) in the eigen directory.
  4. Compile (and link) fndmodlh.for (using the batch file fndlink.bat) in the fndmodlh directory.
  5. Compile (and link) psgraph1.for with psgraph3.for (using the batch file glink1.bat) in the psgraph directory.
  6. Following the instructions in raytracing/samples/readme, run all of the sample cases.
  7. Following the examples in the samples directory, make your own directories with data input files and make the calculations you want.

Directories for the ray tracing program and supporting programs

pcharpo This is the latest version our acoustic ray tracing program for the ocean or atmosphere. It can be compiled in FORTRAN 90. Since January 2006, it can also calculate acoustic-gravity waves if the appropriate dispersion-relation model is used. FORTRAN files: s0.for s1.for s2.for s3.for s4.for conblk2.for (= conblk.for) atmoshd.for ocnhd.for
models This directory includes all of the dispersion-relation models and all of the sound-speed, wind/current, temperature, and gravity models for both the atmosphere or ocean. FORTRAN 90.
fndmodlh Program to search the input data file to determine which models are needed for temperature, etc., and to make a loader file that lists the needed compiled model files. We have batch files that automatically run fndmodlh, and load and run pcharpo. FORTRAN files: fndmodlh.for
eigen Our eigenray program. It calculates the rays that connect a specified transmitter and receiver. It can also plot ground range versus elevation angle of transmission and ground range versus pulse travel time. FORTRAN 90. FORTRAN files: cea.for eigen.for eigen2.for eigen3.for eigen4.for conblk4.for
conplt Plots contours in any vertical or horizontal plane of temperature, sound speed, wind/current velocity component, speed, or direction. FORTRAN 90. FORTRAN files: conplt.for contour1.for conblk3.for
profile Plots vertical profiles of temperature, sound speed, etc. at a specified longitude and latitude. FORTRAN 90. FORTRAN files: profile.for pcprof1.for pcprof2.for pcprof3.for pcprof4.for conblk.for (= conblk1.for)
g This directory includes general-purpose routines that are used by the ray tracing program, the contour-plotting program, the profile-plotting program, and the eigenray program. FORTRAN files: s5.for s6.for s7.for s8.for
psgraph Converts plotting commands output by our other programs to PostScript. FORTRAN 90. FORTRAN files: psgraph1.for psgraph3.for
samples Input and output data files for the sample cases (Look at the samples/readme file for instructions.)

In addition, there are some earlier versions of some of our programs (some for other computers, some for other operating systems, and some for calculating the propagation of radio waves in the ionosphere).

Most of the programs were written by us, but some are modifications of our programs by others. We have never tried to run the programs that are modifications of our programs by others. We furnish them as is, and cannot answer questions about them.

pcharpo.ms models.ms eigen.ms conplt.ms profile.ms psgraph.ms Directories with the extension "ms" are Microsoft FORTRAN versions. NOTICE. If you do not find a model that is mentioned in the HARPA manual or in the HARPO manual in the /models/ directory, look in the /models.ms/ directory. We may not yet have converted the older Microsoft version to a more modern version.
pcharpo.orig This is the original Microsoft FORTRAN version of pcharpo that was published in 1990.
atmospheric_models Some additional atmospheric models. ca. 1980.
harpa.103 Our original acoustic ray tracing program for the atmosphere. Spherical earth. FORTRAN IV for a mainframe computer, 1982.
harpa Acoustic ray tracing program for the atmosphere with irregular terrain. FORTRAN IV for a mainframe computer, 1986.
pcharpa.88 A pc version of HARPA, modification by others, 1988.
harpo Acoustic ray tracing program for the ocean with irregular bottom and an upper surface. FORTRAN IV for a mainframe computer, 1986.
nharpo Revised version of harpo. ca 1986.
3d-iono.75 Our ionospheric radiowave ray tracing program. 1975
pcion-87.apr A pc version of our ionospheric radiowave ray tracing program, modification by others. April 1987
pcion-87.may A pc version of our ionospheric radiowave ray tracing program, modification by others. May 1987
pcion-92 A pc version of our ionospheric radiowave ray tracing program, modification by others. 1992
in_complex_space Ray tracing program in complex space for radio waves in the ionospheric D region for reflection from a conducting layer, in Cartesian coordinates. 1968.
ionospheric_models Additional ionospheric models for ionospheric radiowave ray tracing. ca 1970.

Notice about numerical integration errors and computer word length

The ray tracing program uses an Adams-Moulton predictor-corrector numerical integration method with a Runge-Kutta starter. The user has a choice between testing on relative or absolute single-step error. We have found that relative single-step error works better than absolute single-step error, so we always use that choice, and we recommend that the user make that choice also. With that choice, the following discussion is relevant.

The user has several choices regarding the parameters that control numerical integration errors. Among these, the most important are relative single-step error (W42) and minimum step length (W46). As numerical integration proceeds, the numerical integration subroutine adjusts the step length to guarantee that the actual single-step error for each variable being integrated is less that the requested single-step error (W42). The subroutine keeps decreasing the step length until the required accuracy is achieved. On the other hand, if the actual single-step error is much greater (about a factor of 50, which can be adjusted, but we have never adjusted this) than that required, then the step length is increased.

Problems can arise, however, because of the finite word-length used by the computer to store the variables being integrated. In the most extreme case, if the amount being added to a variable is smaller than the least-significant bit in the computer word, then nothing would actually be added to that variable. But before that extreme case, accuracy could be degraded because of the finite word-length, even when using double precision calculations.

To deal with that problem, we provide a minimum step size, which can be set by the user in the input data file. This is W46, that is, W(46) in the W-array. Clearly, if the minimum step size, W46, is set too small, it will not help deal with the problem. Although, this problem occurs for all of the variables being integrated, it is the height variable (actually the radial variable in the Earth-centered spherical polar coordinate system used by the program) that is most critical because it includes the Earth radius (6370 km).

To help guide the user, we have derived a formula for the minimum value of the product of relative single-step error (W42) and minimum step size (W46). This is:

The product of W42 (maximum allowable single step error) and W46 (minimum integration step size) has to be greater than 10E4 (= approximate Earth radius) / 10E15 (15-digit double precision) This equals 10E-11 km for double precision in FORTRAN 90.

How to find out what the integration error is

When the ray tracing program integrates Hamilton's equations to calculate a raypath, it integrates six equations, three for the position of the ray, and three for the three components of the wavenumber (a vector pointing in the wave-normal direction whose magnitude equals 2 pi divided by the wavelength). The latter three equations are redundant in the sense that given two of the wavenumber vector components, the third can be found from the dispersion relation. However, we integrate all three equations because a comparison with the magnitude of the wavenumber calculated from the dispersion relation gives an estimate of actual numerical integration errors. We calculate this as a relative error, so that it can be compared with the requested single-step error (W42 described above). The program prints out this value in the left-hand column of the standard printout (file DOUTP), and also in the left-hand column of the screen output.

This value should usually be about the same as W42, but when it is not, it usually indicates an inconsistency between the formula for the gradient of some quantity (such as temperature) and the spatial variation of that quantity as specified by the corresponding model. We use this comparison to debug new models.

Sometimes this printed estimate of the error may itself be incorrect if the calculation of the wavenumber from the dispersion relation is wrong. Although this does not happen with any of the acoustic wave dispersion relations, nor with the gravity wave dispersion relation without wind, it sometimes happens for gravity waves with wind. That is, there are cases (discussed below) when the roots of the quartic equation are misidentified so that the wrong root of the quartic equation is chosen, leading to an incorrect value for the wavenumber and therefore to an incorrect estimate of the integration error. When looking at the printout in the DOUTP file, this can usually be recognized because the error may suddenly jump to a very large value, and then later return to a small value.

To deal with this problem, we have generalized the use of the input value of W73, which controls diagnostic printing. For W73=0, there is no diagnostic printing. For W73=1, there is diagnostic printing as described in the HARPO manual. For W73=2 or W73=3, there is additional printing to attempt to print the correct error even when the wrong root of the quartic equation has been chosen.

The usual problem is that sometimes (usually for large wind speed) the roots of the quartic are misidentified, so that the ordinary-ray wavenumber is identified as a small-extraordinary-ray wavenumber, and vice versa. Thus, the wavenumber for the "extraordinary ray" (described below) is chosen when the wavenumber for the "ordinary ray" (also described below) was wanted. If that happens, the program will print out both values (when W73=2 or W73=3), labelling the correct error as "X wave", indicating that the root that had been incorrectly indentified as the "extraordinary" actually has the wavenumber for the "ordinary ray". The incorrect error is also printed.

When W73=3, there is additional diagnostic printing to the screen whenever the program thinks the wrong root has been chosen, or when W73=4 whenever ordinary printing to the standard output file (DOUTP) occurs. In this case, 14 lines are printed to the screen. (This printing to the screen can be redirected to a file to be analyzed later.) What is printed in these 14 lines is described here.

The first column of these 14 lines gives the following 14 quantities:
(1) the relative error, that is, the square of the wavenumber from numerical integration divided by the square of the wavenumber from the dispersion relation minus 1.
(2) the value of the Hamiltonian (which should be zero if the dispersion relation is satisfied) divided by the largest term in the formula for the Hamiltonian.
(3-7) the coefficients of the quartic for (Ck), starting with the coefficient of the constant term, and ending with the coefficient of the fourth power.
(8) the component of the wind in the direction of the wave normal divided by the speed of sound, that whole quantity, squared.
(9-14) the number of the requested ray (0 for the ordinary ray, -1 for the extraordinary ray, -2 for the large extraordinary ray, and 1, 2, 3, or 4 to choose roots 1, 2, 3, or 4).
A description of what is not in the first column is given below.

The first row gives the reason for printing, the height of the ray in km, the ground range of the ray in km, and the values of estimated root numbers for the two ordinary rays (n3 and n4), the large extraordinary ray (nbig), and the small extraordinary ray (nsmall). The numbers 1, 2, 3, and 4 refer to the four columns of complex numbers described below.

The next three rows of eight columns give the complex solutions of the quartic equation. The four pairs of columns correspond to the four complex roots, numbered 1, 2, 3, and 4. The first of the three rows gives the roots actually used. The second row gives approximate roots for small wind. The third row gives the roots from the quartic equation solver (which is not very accurate when the wind is too small).

The next three rows of eight columns give the square of the complex wavenumber from the dispersion relation corresponding to the four complex roots, numbered 1, 2, 3, and 4. The first of the three rows gives the square of the wavenumbers actually used. The second row gives approximate square of the wavenumbers for small wind. The third row gives the square of the wavenumber from the quartic equation solver (which is not very accurate when the wind is too small).

The next three rows of eight columns give the relative error in the solutions to the quartic equation for each of the four complex roots, numbered 1, 2, 3, and 4. Specifically, printed are the values of the quartic for the four roots divided by the maximum term in the quartic. The first of the three rows gives the values of the quartic for the roots actually used. The second row gives the values of the quartic for the approximate roots for small wind. The third row gives the values of the quartic for the roots from the quartic equation solver (which is not very accurate when the wind is too small).

The next three rows of eight columns give the relative error (square of the wavenumber from numerical integration divided by the square of the wavenumber from the dispersion relation minus 1) for each of the four complex solutions of the quartic equation. The four pairs of columns correspond to the four complex roots, numbered 1, 2, 3, and 4. The first of the three rows gives relative error for the roots actually used. The second row gives relative error for the approximate roots for small wind. The third row gives relative error for the roots from the quartic equation solver (which is not very accurate when the wind is too small).

The 14th row gives the square of the wavenumber from numerical integration and the real part of the square of the wavenumber from the dispersion relation.

The gravity-wave dispersion relation with wind

The gravity-wave dispersion relation with wind is a quartic equation. The significance of that is that choosing the frequency and choosing the wave-normal direction is not enough to specify the magnitude of the wavenumber for a gravity wave with wind. There are four possibilities, corresponding to the roots of a quartic equation.

As the magnitude of the wind approaches zero, two of these roots correspond to the wavenumbers for the gravity wave without wind, and the other two correspond to gravity waves whose intrinsic frequency (the frequency as seen by an observer moving with the wind) is very close to the Brunt-Vaisala frequency times the cosine of the angle of the wavenormal direction relative to the horizontal. The wavenumbers for these latter two roots approach infinity as the wind approaches zero, one of the roots approaching infinity much faster than the other. (That is, the wavelengths of these latter two roots approach zero as the wind approaches zero.)

We refer to the former two roots as the "ordinary ray." For small wind, these correspond to a wave whose wavenormal points in the direction we wish to propagate, and to another wave whose wavenormal points in the opposite direction. Of the latter two roots, we refer to the root that approaches infinity faster as the wind approaches zero as the "large extraordinary ray," and the other root as the "small extraordinary ray." Notice that these different kinds of rays are analogous to ordinary and extraordinary waves for radiowave propagation in the ionosphere.

Since these different kinds of waves have different properties, it is important to be able to distinguish among them, not only so that the raypath calculated is actually the raypath wanted (because this choice determines the initial conditions for a raypath calculation), but also so that estimated errors are correctly calculated, as discussed above.

We have also found a case where a ray that started out as an "ordinary ray" transformed into a "small extraordinary ray." As far as we could tell, this was an actual transformation rather than a misidentification of the roots of the quartic. As the ray was propagating, roots of the quartic for the "ordinary ray" and the "small extraordinary ray" became closer. When they became equal (within a percent), the ray tracing program began calculating the "small extraordinary ray" rather than the "ordinary ray." However, there was no discontinuity. Everything was smooth. The relative error increased only slightly at the first integration step for the "small extraordinary ray," then recovered by the next integration step. The process seems to be analogous to waveguide mode coupling. Here is an example.

Documentation/hgrav18_coupling_example. In this example, we have an isothermal atmosphere of temperature 224 K. We have a wind blowing from West to East with a constant wind shear of 0.1 m/s/km, with zero wind at the ground. The plot projection is on a vertical East-West plane with West on the left.

The first frame shows downwind propagation of the "ordinary ray" for elevation angles from -60 degrees through +60 degrees. Plus and minus 70 and 80 degrees elevation angles of transmission are evanescent.

The second frame shows upwind propagation of the "ordinary ray" for elevation angles of transmission of -60 degrees through +60 degrees. Plus and minus 70 and 80 degrees elevation angles of transmission are evanescent. The + and - 60 degree elevation angle of transmission (the ones with the big loops) are actually examples of coupling between the "ordinary ray" and the "small extraordinary ray." The coupling occurs at a height of about 3 and a half km. We have the "ordinary ray" below that height and the "small extraordinary ray" above that height, so the coupling goes both ways. The wind shear gives a clockwise vorticity, which is the direction of ray travel in the large loops.

The third frame shows downwind propagation of the "small extraordinary ray" for elevation angles of plus and minus 70 and 80 degrees. Other angles of transmission are evanescent. The rays are very close to the ground.

The fourth frame shows upwind propagation of the "small extraordinary ray" for elevation angles of plus and minus 60 degrees. Plus and minus 70 and 80 degrees elevation angles of transmission are evanescent. Other angles of transmission do not get very high above the ground.

Below is the standard printout for these rays.
Documentation/Dgrav18_coupling_example. The first non-evanescent ray is for downwind propagation of the "ordinary ray" for an elevation angle of transmission of -60 degrees. This has an upward ray direction. There is a misidentification of the ray from about 1.7 km to 2.6 km, as can be seen. This is not coupling. The actual "small wavenumber ray" is evanescent at this height. The +60 degree elevation angle of transmission behaves similarly.

The first non-evanescent ray for upwind propagation of the "ordinary ray" is for -60 degrees elevation angle of transmission. At a height of about 3.7 km, it couples to the "small extraordinary ray," where it remains until returning to about that height, when it couples back to the "ordinary ray," reflects from the ground, again couples to the "small extraordinary ray" at a height of about 3.7 km, again goes up to about 333 km, comes back down, and couples back to the "ordinary ray" at a height of about 3 and a half km before hitting the ground. The behavior of the +60 elevation angle of transmission case is similar, except that it reflects off the ground first, and then makes only one big loop.

For downwind propagation of the "small extraordinary ray," it can be seen that only elevation angles of transmission of plus and minus 70 and 80 degrees propagate (although very close to the ground), while other directions are evanescent.

For upwind propagation of the "small extraordinary ray," the first non-evanescent ray is for -60 degrees elevation angle of transmission. It makes two large clockwise loops. The last non-evanescent ray is for +60 degrees elevation angle of transmission. It makes one large clockwise loop. The other rays do not get very high as can be seen.

The diagnostic printing to the screen was redirected to the following file.
Documentation/testing18_coupling_example. There are seven columns. The first column gives the relative error between the integrated wavenumber and the wavenumber calculated from the dispersion relation. In most cases here, large errors in this column are due to a misidentification of the quartic roots. In such cases, the program also prints out the root with the smallest error, and prints out the identification given to that root by the dispersion relation subroutine. The second column gives that identification. The third column gives the height in km; the fourth column gives the ground range in km. The fifth column gives the value of the root that has been identified as the "ordinary ray." The root equals the square of the wavenumber times the square of the sound speed. The sixth column gives the value of the root that has been identified as the "small extraordinary ray." The sixth column gives the elevation angle of transmission.

The first propagating ray is for downwind propagation of the "ordinary ray" for an elevation angle of transmission of -60 degrees. As can be seen, the roots were misidentified from a height of about 1.7 km to 2.6 km. The "small extraordinary ray" is evanescent (negative square of wavenumber) at those heights.

The first propagating ray for upwind propagation of the "ordinary ray" is for an elevation angle of transmission of -60 degrees. As can be seen, coupling from the "ordinary ray" to the "small extraordinary ray" occurs at a height of about 3.7 km, where the square of the "ordinary ray" wavenumber and the square of the "small extraordinary ray" wavenumber are about equal.

The input data file for these ray calculations is at
Documentation/grav18_coupling_example.

Dispersion relation for magneto-acoustic-gravity waves.

We derived the dispersion relation for magneto-acoustic-gravity waves. A reprint of the paper published in JASTP is at http://cires.colorado.edu/~mjones/raytracing/Publish/MagnetoAcousticGravityWaveFigures/1-s2.0-S1364682616303406-main.pdf . We hope to eventually use that dispersion relation to make a dispersion-relation model for the ray tracing program.

Algorithms we use to find intersections of the ray with the terrain and with the receiver surface

The algorithms we use are explained in our report Documentation/NOAA_WPL98.pdf.

Guidelines we have tried to follow in developing our ray tracing programs

  1. The program should be easy to modify for specific purposes, not necessarily so general to begin with. For this, we need a good division of labor between subprograms so that it is necessary to change only one or two subprograms for specific purposes.
  2. The input to the W array should be such that input data are ready to use directly, without having to convert or change units every time the variable is used. Thus, SUBPROGRAM READW converts angles to radians as soon as they are read in, for example.
  3. We should try not to calculate something many times, but instead, we should put the variable in labeled common and calculate it outside of subprograms and loops.
  4. We should use actual variable names in blank common even though the variable names may be equivalenced to an element in an array so that there is a closer connection to the variables used in the physical situation.
  5. All subprograms that read in data should print it back out with headings.

Updates to the ray tracing programs

  1. We have combined some subprograms that are now shared among the various programs. These are in the subdirectory /raytracing/g/.
  2. We have updated the eigenray program to allow ends of rays to be used as virtual transmitters to calculate Fresnel zones. The eigenray program produces a file named doutp2, that can be used as input to the ray tracing program.
  3. We have fixed the problem of raypaths occasionally entering the plotting region at the wrong place. (19 April 2007)
  4. Log base 10 of density in kg per cubic meter is now output to the transmitter and receiver raysets in columns 81-90 with an implied decimal point between columns 85 and 86. (20 April 2007)
  5. Added signal strength calculations (including spreading loss) to the eigenray program (eigen). The spreading loss calculations are valid for an atmosphere that is azimuthally symmetric about the transmitter. That condition is satisfied by an atmosphere that varies only with height. The output is to a file named doutp3. (1-7 May 2007)
  6. The Eigenray program will now allow raysets generated by the ray tracing program in which frequency is stepped. (4-8 May 2007)
  7. Added output to doutp4 and doutp5 in the eigenray program. (23 Nov 2007)
  8. With the help of Larry Baker of USGS, we have fixed some common block length discrepancies that will allow our programs to run on a larger variety of systems. (15 February - 19 March 2008)
  9. With the help of Larry Baker of USGS, we have corrected various problems, that while not affecting the calculations we have made, will allow our programs to run on a larger variety of systems. (15 February - 19 March 2008)
  10. We now use a different routine for date and time (courtesy of Larry Baker of USGS), (5 March 2008)
  11. We have updated the calculation of machine constants using I1mach.f90, R1mach.f90, and D1mach.f90 (26 August 2008) (D1mach is used in vdraft2, s5, contour1, pcprof1, and pcprof2.)
  12. With the help of Larry Baker of USGS, we have fixed some type discrepancy problems and array dimensioning problems (10 September 2008)
  13. Fixed common block length discrepancies found by Larry Baker at USGS that will allow the ray tracing program and the other programs to run on a larger variety of systems. (September 2008)
  14. Fixed other incompatibility problems found by Larry Baker at USGS that will allow the ray tracing program and the other programs to run on a larger variety of systems. (September 2008)
  15. Added units conversion for a central Earth angle in nautical miles. (19 Mar 2009)
  16. To improve portability of the program we replaced REAL*8 by DOUBLE PRECISION (30 Apr - 4 May 2009).
  17. We added plotting of a circle to projection on the ground showing the core of the vortex whenever the VVORTX3 wind model is used - (16-29 June 2009).
  18. We added optional labels to the projection on the ground giving some of the parameters of the vortex whenever the VVORTX3 wind model is used. Documentation is in SUBROUTINE PLOTNOTE in file s3.for - (18-29 June 2009).
  19. We added the possibility of increasing the size of the above labels and increasing or decreasing the size of other labels as well. The size variable can have the values 0, 1, 2, or 3, and corresponds to multiplying the default font size by 0.75, 1.0, 1.5, or 2.0. This involved changes in S3.for, S6.for, S8.for, contour1.for, pcprof3.for, and PSGRAPH1.FOR (23 June 2009)
  20. The eigenray program (eigen) now outputs time, phase, absorption, and Doppler shift to the file doutp2. (13 August 2009)
  21. The eigenray program (eigen) now interpolates phase, absorption, and Doppler shift in addition to pulse travel time. (13 August 2009)
  22. To make calculating Fresnel zones easier using the eigenray program, W(453) now takes the desired receiver range as a central earth angle in radians, degrees, or km. W(459) is now an input to the eigenray program as a phase to subtract from interpolated phases. For Fresnel zone calculations, this would normally be the phase of the direct raypath. (14 August 2009)
  23. Added more documentation to the output file doutp3 in the eigenray program. (29 October 2009)
  24. Corrected the degree tick mark values in cross range for plot projections on the ground in SUBROUTNE PLTANH in file s6. (3 Nov 2009)
  25. Used D1mach instead of rmach for machine constants in VDRAFT2 (6 Nov 2009)
  26. Our contour program (conplt) can now label contours on contour plots. It is not automatic, but it can be done using the input data, which allows us to produce publication-quality contour plots without having to do manual post processing. Some examples are at Documentation/ctilt0.pdf. The input data file for producing these contour plots is at Documentation/tilt0. (18 March 2010)

    This labeling uses a new subroutine PLTNOTE, which is in s5.for. PLTNOTE can be used by other programs also for putting labels on plots.

  27. Made sure the entry points were consistent in the gravity-wave dispersion relation model and in the pressure, density, and gravity models. (22-29 March 2010)
  28. Added a new dispersion-relation model for gravity waves with wind, no losses. (29 April 2010) See new models for details. .
  29. Fixed a bug that labeled the height axis in plots as being in km even when the heights were actually in meters. The correction was in s6.for (2 November 2010)
  30. Generalized the formulas for vertical wavenumber in all of the dispersion relation models to handle the case where the reference height (receiver height) is in an evanescent region. These formulas are used for profile plotting. (7-13 April 2011)
  31. Generalized the formulas for wavelength and phase speed in the profile program (file pcprof3.for) to handle the case where the reference height (receiver height) is in an evanescent region. (7-13 April 2011)
  32. Added the ability to put a plus sign on a raypath at the place where the pulse travel time reaches a specified value. The value wanted is specified in seconds in W78. (Files s1.for, s3.for, and s6.for required modification.) (13 April 2011)
  33. In the profile program, for terrain profiles, it is now possible to choose the minimum height to be w(88), the maximum height to be w(89), and the distance between vertical tic marks to be w(96) instead of getting auto scaling. The left latitude w(83) and right latitude w(85) will determine the left and right edges of the plot, but not necessarily in that order. Because of auto scaling for the horizontal axis, colatitude will increase from left to right. If the left latitude is always chosen to be larger than the right latitude, then the profile plot can be made to match a raypath plot. (22 September 2011)
  34. Modified the diagnostic printing in SUBROUTINES TRACE and PRINTR for the gravity-wave dispersion-relation model with winds and no losses. (Files s1.for and s2.for) (2 March 2012)
  35. Modified the dispersion-relation model for gravity waves with wind, no losses. (2-9 March 2012) See new models for details. .
  36. Modified the dispersion-relation model for gravity waves with wind, no losses to improve identification of the roots of the quartic equation that determines the magnitude of the wavenumber. (28-29 May 2013) See new models for details. .
  37. Adding Coriolis terms to the two gravity-wave dispersion relation models. This is being done in steps:
    1. Added Coriolis terms to the /GV/ common block. This affects S6, GWWNL, GNWNL, NGRAV, GCONST, and NPGRAV. (27 June 2016)
    2. Changed calling sequence from (gravity) to (gravity,coriolis) for gravity models. (27 June 2016)
    3. Added calculations of the Earth's angular velocity to gconst. (28 June 2016)
    4. Added the main Coriolis term to the Hamiltonian in the gravity wave dispersion relation without wind GNWNL. (29 June 2016)
    5. Added the main Coriolis term to the calculation of kay2 and kz2a in the gravity wave dispersion relation without wind GNWNL. (30 June 2016)
    6. Added the four more Coriolis terms to the Hamiltonian in the gravity wave dispersion relation without wind GNWNL. (22 July 2016)
    7. Added the four more Coriolis terms to the calculation of kay2 and kz2a in the gravity wave dispersion relation without wind GNWNL. (22 July 2016)
    8. Have tested, debugged, and corrected the above four terms. The main difficulty was getting the correct formulas for changing the normal component of the wave-normal vector at reflection from the terrain. However, the present formulas with Coriolis terms work only for a horizontal surface. (9 August 2016)
    9. Added the 6th (and last) Coriolis term to the gravity wave dispersion relation without wind GNWNL. (14 September 2016)
    10. Added all 6 Coriolis terms to the gravity wave dispersion relation with wind GWWNL. (21 July 2017)
      This involved several steps to complete this addition (2 August 2017):
      • Checking the integration error, and correcting the formulas in the added Coriolis terms until the error is no greater than the requested error (w42).
      • Checking the error calculation to make sure it is correct for the added Coriolis terms.
      • Checking the formulas for reflection from the terrain, and making sure they are correct for the added Coriolis terms.
      • Checking the approximations for ordinary and extraordinary rays to make sure they are correct for the added Coriolis terms.
  38. Added labels to raypath plots using PLTNOTE. This involved an addition to S3.for. The labels go into the data for the description of the wind model and the description of the wind perturbation model. There are a maximum of 16 labels. The x and y positions for the labels go into W101-w116 and w126-141. Set w100=-6 and w125=-6 to use. These labels can be used only for ray tracing models that do not have wind. For ray tracing models that have wind, something similar will have to be added. An example of usage is in the file
    (http://cires.colorado.edu/~mjones/raytracing/antarctica/grav11
    . (14 Dec 2016)
  39. Improved the calculation of wavelength and vertical wavenumber in the profile program to handle the extra complications of the gravity-wave dispersion relation with wind, especially for cases of extraneous waves. This involved adding a new common block /RINK4/. The calculation for profiles of phase velecity and phase integral will also have to be modified in the same way. (25 June - 9 July 2018)

New models in subdirectory /raytracing/models/

  1. A wind profile whose magnitude and direction can approximate measurements arbitrarily well
    (http://cires.colorado.edu/~mjones/raytracing/models/vtanh.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/Vtanh.pdf
    ) (18 June 2004)
  2. An updraft/downdraft wind perturbation model with inflow/outflow at the ground (
    (http://cires.colorado.edu/~mjones/raytracing/models/vdraft.for
    , described by http://cires.colorado.edu/~mjones/raytracing/models/Vdraft.pdf ). (June 2005)
  3. A new dispersion relation for acoustic-gravity waves, no wind, no losses
    (http://cires.colorado.edu/~mjones/raytracing/models/gnwnl.for)
    . Using this dispersion-relation model requires temperature models that furnish second derivatives of the temperature. To do this, we added additional entry points in the temperature models to furnish the temperature plus first and second derivatives. The original entry points still furnish only temperature plus first derivatives for use with acoustic waves. (January-March 2006)
  4. A pressure model based on a varying scale height that depends on temperature ( http://cires.colorado.edu/~mjones/raytracing/models/pscaleh.for ). (January-March 2006)
  5. A constant gravity model ( http://cires.colorado.edu/~mjones/raytracing/models/gconst.for ). (January-March 2006)
  6. A generic gravity model ( http://cires.colorado.edu/~mjones/raytracing/models/ngrav.for ). (January-March 2006)
  7. A generic gravity perturbation model ( http://cires.colorado.edu/~mjones/raytracing/models/npgrav.for ). (January-March 2006)
  8. A simple density model based on hydrostatic equilibrium and a given temperature model
    (http://cires.colorado.edu/~mjones/raytracing/models/rhosmp.for
    ). (January-March 2006)
  9. A generic density model ( http://cires.colorado.edu/~mjones/raytracing/models/nrho.for ). (January-March 2006)
  10. A generic density perturbation model ( http://cires.colorado.edu/~mjones/raytracing/models/nprho.for ). (January-March 2006)
  11. A temperature perturbation model based on the stream function for the updraft/downdraft wind perturbation model with inflow/outflow at the ground
    (http://cires.colorado.edu/~mjones/raytracing/models/tdraft.for
    , described by http://cires.colorado.edu/~mjones/raytracing/models/Tdraft.pdf ). (April 2006)
  12. An updraft/downdraft wind perturbation model with inflow/outflow at the ground and outflow/inflow at a specified height
    (http://cires.colorado.edu/~mjones/raytracing/models/vdraft2.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/Vdraft2.pdf
    ). (August 2006)
  13. A temperature perturbation model based on the stream function for the updraft/downdraft wind perturbation model with inflow/outflow at the ground and outflow/inflow at a specified height
    (http://cires.colorado.edu/~mjones/raytracing/models/tdraft2.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/Tdraft2.pdf
    ). (August 2006, revised November 2006)
  14. A new loss model based on the formulas in the paper by Bass and Hetzer (Inframatics No. 15, September 2006, pp. 1-12) (
    (http://cires.colorado.edu/~mjones/raytracing/models/sbloss.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/SBloss.pdf
    ). (March 2007)
  15. The constant molecular-weight model ( http://cires.colorado.edu/~mjones/raytracing/models/mwconst.for ) is generalized to include the amount of N2, O2, N, O, H2O, and O3. The generalization is described in http://cires.colorado.edu/~mjones/raytracing/models/MWconst.pdf . (March 2007)
  16. A new temperature perturbation model to approximate the temperature distribution inside of a hurricane
    (http://cires.colorado.edu/~mjones/raytracing/models/Tplume.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/Tplume.pdf
    ) (10 April 2009)
  17. A new wind/current perturbation model of a vortex to represent a hurricane
    (http://cires.colorado.edu/~mjones/raytracing/models/vvortx4.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/Vvortx4.pdf
    ) This model should be used instead of VVORTX3 for a vortex anyplace other than the equator. VVORTX4 gives the same vortex no matter where it is put, but VVORTX3 gives a different vortex when moved off of the equator. For example, when using VVORTX3, we get the following raypaths for a vortex at 40 degrees north latitude. hurrican/Hurricane at 40 deg N using vvortx3.pdf . On the other hand, when using VVORTX4, we get hurrican/Hurricane at 40 deg N using vvortx4.pdf , which is the same as either VVORTX3 or VVORTX4 would get if the vortex were on the equator. (1 July 2009)
  18. An updraft/downdraft wind perturbation model with inflow/outflow at the ground and outflow/inflow at a specified height with a vortex ring at the top
    (http://cires.colorado.edu/~mjones/raytracing/models/vdraft3.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/Vdraft3.pdf
    ). (December 2009)
  19. Added a new dispersion-relation model for gravity waves with wind, no losses.
    http://cires.colorado.edu/~mjones/raytracing/models/gwwnl.for
    , (April-June 2010). An example of the raypaths is given at
    http://cires.colorado.edu/~mjones/raytracing/Documentation/hgrav17sample.pdf
    . Another example of the raypaths is given at
    http://cires.colorado.edu/~mjones/raytracing/Documentation/hgrav17test3Jun2013.pdf
    .

    When winds are added to gravity waves, it is necessary to solve a quartic equation to calculate the wavenumber (two pi divided by the wavelength) from knowing the frequency of the wave and the direction of propagation. This means there is the possibility for four different kinds of gravity waves propagating in the same direction, similar to ordinary and extraordinary waves for radiowaves in the ionosphere. However, some of these solutions to the quartic may be extraneous. For example a negative value for the wavenumber would be extraneous because we have already defined the wavenumber as a positive magnitude. (Actually, a negative nalue for the magnitude of the wavenumber in this case corresponds to a wave propagating in the opposite direction from the wavenormal direction, but it is still extraneous.)

    Two of these waves correspond to the usual gravity waves in the absence of wind. In analogy with radio waves, we refer to these as ``ordinary'' waves. One of these has a negative value, and is therefore extraneous.

    The other two waves are not present without wind. Their wavenumbers approach infinity as the wind speed approaches zero. Again, in analogy to electromagnetic waves, we refer to these waves as ``extraordinary'' waves. For small wind speed, the wavenumber for one of these waves is much larger than the wavenumber for the other. We refer to the former wave as a ``large extraordinary'' wave, and we refer to the latter as a ``small extraordinary'' wave, or sometimes as simply an ``extraordinary'' wave. These waves correspond to waves whose intrinsic frequency (the frequency as seen by an observer moving with the wind) is approximately equal to the Brunt-Vaisala frequency times the cosine of the elevation angle of the wave normal.

    The model now more accurately identifies the four roots of the quartic equation, but still makes some errors when the wind speed is too large. (28-29 May 2013). See below for details. .

    In the example raypaths shown here, the first frame shows examples of ordinary rays. The second frame shows an example of a small extraordinary ray.

    The input data is at
    http://cires.colorado.edu/~mjones/raytracing/Documentation/grav17
    . The printout showing the details are at
    http://cires.colorado.edu/~mjones/raytracing/Documentation/Dgrav17
    .

    The input data for the second example is at
    http://cires.colorado.edu/~mjones/raytracing/Documentation/grav17test3Jun2013
    . The printout showing the details are at
    http://cires.colorado.edu/~mjones/raytracing/Documentation/Dgrav17test3Jun2013
    .

  20. The temperature perturbation model Tplume is now generalized to allow for two plumes.
    (http://cires.colorado.edu/~mjones/raytracing/models/Tplume.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/Tplume.pdf
    ) (9 September 2010)
  21. The temperature perturbation model Tplume is now generalized to allow for more abrupt edges for the plumes. In addition to a Gaussian fall-off with distance from the center of the plume, we also provide for an exponential dependence on the 4th, 6th, and 8th power of the distance from the center of the plume.
    (http://cires.colorado.edu/~mjones/raytracing/models/Tplume.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/Tplume.pdf
    ) (29 October 2010)
  22. Added a new temperature perturbation model for a temperature inversion that varies as a harmonic function of latitude, as though it were modulated by a gravity wave.
    (http://cires.colorado.edu/~mjones/raytracing/models/Twave.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/Twave.pdf
    ) (27 December 2010)
  23. The temperature model Tlinear is now generalized to be used for both gravity waves and acoustic waves.
    (http://cires.colorado.edu/~mjones/raytracing/models/Tlinear.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/Tlinear.pdf
    ) (16 March 2011)
  24. The terrain (bottom) model Gtanh has been converted from MicroSoft FORTRAN to FORTRAN 90.
    (http://cires.colorado.edu/~mjones/raytracing/models/gtanh.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/Gtanh.pdf
    ) (19 August 2011)
  25. We have a new terrain (bottom) model Gtanh2, which is equivalent to Gtanh, except that it uses the general subrouting FTANH2 to calculate the profiles.
    (http://cires.colorado.edu/~mjones/raytracing/models/gtanh2.for
    , described by
    http://cires.colorado.edu/~mjones/raytracing/models/Gtanh2.pdf
    ) (19 August 2011)
  26. Modified the dispersion-relation model for gravity waves with wind, no losses http://cires.colorado.edu/~mjones/raytracing/models/gwwnl.for (2-9 March 2012) to allow choosing any of the four roots of the quartic (by using W2 = 1, 2, 3, or 4 for roots 1, 2, 3, and 4) for the initial wavenumber at the transmitter. We also now use W2 = 0 for the ordinary wave, W2 = -1 for the small extraordinary wave, and W2 = -2 for the large extraordinary wave.
  27. Modified the dispersion-relation model for gravity waves with wind, no losses http://cires.colorado.edu/~mjones/raytracing/models/gwwnl.for (28-29 May 2013) to more accurately identify the roots of the quartic equation to distinguish among ordinary ray, small extraordinary ray, and large extraordinary ray.
  28. A full list of the models is given here.

Profile and contour plotting

  1. Sound speed (w90 = 1)
  2. Wind speed (w90 = 2)
  3. Vertical component of wind (w90 = 3)
  4. Southward component of wind (w90 = 4)
  5. Eastward component of wind (w90 = 5)
  6. Topography (w90 = 6)
  7. Temperature (w90 = 7)
  8. Wind direction (w90 = 8)
  9. Acoustic cutoff frequency (w90 = 9)
  10. Brunt-Väisälä frequency (w90 = 10)
  11. Square of Brunt-Väisälä frequency (w90 = 11)
  12. Vertical wavenumber (w90 = 12)
  13. wavelength (w90 = 13)
  14. phase speed (w90 = 14)
  15. phase integral (w90 = 15), profile only, no contours (to show regions of validity of the WKB approximation)

Improvements still to be made to the ray tracing program

  1. Add Coriolis force to the two gravity-wave dispersion relation models. This will require several steps:
    1. Need to change the common block /GV/ that has the output data from the gravity models to include the magnitude of Omega, its components in the r, theta, and phi directions, and the derivatives of all 4 of those with respect to time, r, theta, and phi. In total, 20 additional elements. (DONE)

      All of the routines that refer to the /GV/ common block will have to be revised. These are S6, the gravity wave dispersion relations GWWNL and GNWNL and the gravity models NGRAV and GCONST and the gravity perturbation model NPGRAV. (DONE)

    2. Need to change the calling sequences from (gravity) to (gravity, coriolis) for the two gravity models NGRAV and GCONST and the gravity perturbation model NPGRAV. (DONE)
    3. All of the routines that call the gravity models will have to be revised. These are the pressure model PscaleH and the density model RHOSMP. (DONE)
    4. the gravity model GCONST (that now calculate the Earth's gravitational field and it's gradient) will have to be revised to calculate in addition the components of the Earth's angular velocity and all components of their gradients and add all of that to the output from these models. (DONE)
    5. Need to add the coriolis terms to the dispersion relation (which is then used as the Hamiltonian) to the gravity-wave dispersion relation model without wind GNWNL. (DONE)
    6. It will also be necessary to calculate the derivatives of those Coriolis terms with respect to r, theta, phi, time, kr, ktheta, and kphi to get refraction and Doppler shift calculations correct for the gravity-wave dispersion relation model without wind GNWNL. (DONE)
    7. Testing and de-bugging the changes may take several months. The first 5 Coriolis terms have been tested and debugged. However, reflection formulas are now valid only for a horizontal surface when Coriolis terms are present. (9 Aug 2016)
    8. Need to add Coriolis term 6 to the gravity-wave dispersion relation model without wind GNWNL. (DONE)
    9. Need to generalize the reflection formulas so that they are correct for arbitrary terrain when Coriolis terms are present.
    10. Need to add the coriolis terms to the dispersion relation (which is then used as the Hamiltonian) to the gravity-wave dispersion relation model with wind GWWNL.
    11. It will also be necessary to calculate the derivatives of those Coriolis terms with respect to r, theta, phi, time, kr, ktheta, and kphi to get refraction and Doppler shift calculations correct for the gravity-wave dispersion relation model GWWNL with wind.
    12. Testing and de-bugging the changes may take several months.
    Details of the dispersion relation including the Coriolis terms are at :
    http://cires.colorado.edu/~mjones/raytracing/Coriolis.pdf .
  2. (2015) Dispersion relation for magneto-acoustic-gravity waves.

    We derived the dispersion relation for magneto-acoustic-gravity waves. A preprint of the paper we plan to submit to JASTP is at http://cires.colorado.edu/~mjones/raytracing/Publish/WavesJASTP/Waves.pdf . We hope to eventually use that dispersion relation to make a dispersion-relation model for the ray tracing program.

  3. The algorithms for reflection from the terrain for gravity waves with wind are valid only for level terrain or horizontal wind. They need to be generalized.
  4. Generalize the spreading-loss calculations in the eigenray program to allow arbitrary horizontal gradients in the atmosphere. This would probably take 1 month working at 21% time.
  5. Add a dispersion relation model for Lamb waves. This would probably take 4 months working at 21% time.
  6. Add the capability of having the transmitter in an evanescent region. There are two possibilities. We could do a theoretical calculation based on the formulas for the WKB approximation (when valid) or using Airy functions (otherwise). That might take a week or less. The other possibility is do ray tracing in complex space. This requires a significant enhancement to the ray tracing program or a special-purpose ray tracing program to do ray tracing in complex space. See the discussion of ray tracing in complex space further down the list.
  7. Some of the model subroutines still need to be updated to be compatible with the upgrading of the ray tracing program for acoustic-gravity waves.
  8. Some of the input data files need to be updated to be compatible with gravity-wave models that require second derivatives of temperature. This would probably take 1 month working at 21% time.
  9. Need to increase the space for printing the numbers on some axis labels. This is from a difference in Microsoft Fortran and Fortran 90. We shall do this as needed.
  10. The wind gradient is discontinuous on the axis of the updraft/downdraft models (vdraft2.for) and (vdraft3.for). It may be a good idea to replace the present updraft/downdraft models with ones that are continuous sometime in the future. This would probably take 2 weeks working at 21% time.
  11. We might add the possibility for atmospheric models to have parameters that depend on height, longitude, or latitude in very general ways, with general subroutines to furnish those parameters to any models that want them. We shall do this as needed.
  12. It might be necessary to add ray tracing in complex space to calculate the propagation of evanescent waves. We would either enhance the ray tracing program to do ray tracing in complex space or develop a special purpose program to do ray tracing in complex space. It would take careful consideration before starting on either of those projects. It might take about a year at 21% time.
  13. We might need to write alternate versions of SUBPROGRAM HAMLTN (that defines the differential equations being integrated) for specific cases. We shall do this if we need to.
  14. We might want to calculate polarization (that is, the amplitudes and phases of perturbation values of temperature, pressure, density, velocity, etc. associated with the wave). We shall do this if we need to.
  15. We might want to integrate focusing. This would probably take about 1 year working at 21% time. It would require adding an additional derivative to all calculations. That is, if we are already taking second derivatives of some quantity, we would need to add third derivatives. This is probably not worth doing.
  16. We might want to simultaneously plot raypaths on a vertical and horizontal plane. This would require a great deal of work. It is probably not worth the effort.
  17. We might want several versions of PRINTR (what and how much to print). We shall do this if we need to.
  18. We might want several versions of TRACE (when and how often to print). We shall do this if we need to.
  19. We might want a simple version of TRACE. We shall do this if we need to.
  20. We might want to alter the PROFILE and CONPLT programs to include conversions between the geographic coordinate system and the computational coordinate system. Possibly. Maybe not.
  21. We might want to include ground reflections in the eigenray program. We shall do this if we need to. It means including ground reflection coefficients that might depend on location and angle of incidence.
  22. In adding together the contributions from several hops to get the total field in the eigenray program, we have to choose a ground range and interpret between raysets. We shall do this if we need to.
  23. In SUBROUTINE TRACE, there is a place where the program divides by D2Z, which could lead to trouble if D2Z is zero, and the program does not test for that possibility. This is in line TRWE1940 on page 273 of the HARPO manual.

Improvements to the documentation of the ray tracing program

Improvements done

  1. We now have a list of all models with their descriptions. It is at models.
  2. We now have a list of all of the routines along with a description of what each one does. It is at Documentation/subroutines.html.
  3. Every routine contains a descriptive comment at its beginning describing what the routine does.

Improvements still to be done

  1. Need to add a .pdf file to describe each of the models. The models that still need .pdf files (including all 6 of the dispersion relation models) are: Those models marked by an asterisk * have the highest priority. Those models beginning with n are do-nothing models.
  2. Still need to add comments to many of the routines to define its major variables, especially its input and output parameters.
  3. Need to make sure the instructions for finding eigenrays are clear.
  4. Need to clearly tell which modules are system dependent (like graphics) and what the user has to do about it.
  5. Include an index.
  6. Include figures from Walter Dudziak's report (used with permission) to show the relationship between the geographical coordinate system and the computational coordinate system.

How to add new models to the ray tracing program

The versatility of this ray tracing program comes mostly from the ability to add new models for the various atmospheric parameters, including temperature, wind, and the dispersion relation. To do that, the easiest way is to take one of the present models that is most like the new one you want, copy it, rename it, and modify it to suit your needs.

However, there are some necessary details that you need to know. First, each model for each type (pressure, temperature, etc.) is assigned a number. If you add a new model, you need to assign a new number, and add that number to the file raytracing/fndmodlh/modltabl.dat .

Another is entry points. These are explained here.

Another consideration is the input data used by a model. Each kind of model (temperature perturbations, for example) is assigned 25 elements in the W array. This assignment is given by Table 4.4 on page 67 in the HARPA manual, or by Table 4.4 on page 66 in the HARPO manual. However, the above tables do not include the models we have added for gravity waves since those reports were published in 1986. These are given in
raytracing/fndmodlh/modltabl.dat
.

In addition, there is allowance for voluminous data (such as profiles) into each model. Each type of model is assigned a common block. The names of these common blocks are also given in raytracing/fndmodlh/modltabl.dat . The structure of these common blocks (using the temperature common block as an example) is given in Table 4.5 on page 69 in the HARPA manual, or (using the sound-speed common block as an example) is given in Table 4.5 on page 68 in the HARPO manual.

Each model also puts its output to a common block. The assignment of common blocks is given in Table 4.3 on page 63 of the HARPA manual, or in Table 4.3 on page 62 in the HARPO manual. In addition, for the new kinds of models, we have the following assignment: Common block /GV/ for gravity, and /DS/ for density.

In addition, there are other considerations when designing a model that will work correctly and efficiently with the ray tracing program.

Each model is implemented by a FORTRAN subroutine that must specify some parameter (such as temperature, sound speed, density, pressure) as a continuous function of height, latitude, and longitude (more specifically, as a continuous function of position in an Earth-centered spherical polar coordinate system). In addition, the subroutine must specify first derivatives (and sometimes second derivatives) as continuous functions of position. In addition, the specification of the derivatives must be mathematically consistent with the spatial variation of the corresponding parameter. Models that fail to follow this rule will cause the ray tracing program to run slowly and give inaccurate results.

Bottom models and upper surface models specify a continuous function that is zero on the surface, positive in the region where the ray is allowed to be, and negative in the region where the ray is not allowed to be. That is, the bottom model function is positive above the bottom. The upper surface model function (for underwater ray tracing) is positive below the upper surface. These models must specify both first and second derivatives as continuous functions of position in Earth-centered spherical polar coordinates. This method of specifying the bottom models allows for overhangs and caves. For atmospheric ray tracing, bottom models would be called terrain models. For underwater ray tracing, they would be called bathymetry models.

Receiver-surface models specify a continuous function that is zero on the receiver surface. These models must specify both first and second derivatives as continuous functions of position in Earth-centered spherical polar coordinates. Notice that a receiver-surface model can be specified as being relative to the terrain (the bottom). The receiver surface model RBOTM is relative to the terrain.