Geometric height calculation discrepancy (ROPP cf GRAS NRT)

[Huw Lewis]

Email from Kjartan:

I think I found a discrepancy between the geopotential to geometric
altitude conversion in  ropp_fm_obs2roprof and the one that was 
used to generate the most recent version GRAS SAF nrt data.

The simplest way to generate the issue is this: 

Read a refractivity profile from the GRAS SAF nrt data (atm* - file)
pass the resulting RO_data structure through 
ropp_fm_roprof2obs(RO_data,obs) and pass it back through 
ropp_fm_obs2roprof(obs,RO_data). 

The RO_data structure now has a vector of geometric altitudes 
(RO_data%Lev2a%alt_refrac) that is inconsistent with the values in
the input file, although all other data are the same. I suspect the
problem derives from a discrepancy in the value used for the 
acceleration of gravity.

Huw's initial response:

You are correct in saying that the input alt_refrac data after passing
through ROPP processing is not identically the same as that in the input
file. You would also notice if you were to do the same for the bending
angle observations that only a limited sub-set of Level1b data from the
input file is preserved on output. Both these characteristics result
from the ropp_fm_obs2roprof defining a new Level2a (or Level1b)
structure to output the (potentially updated) values. I think this is a
valid approach.

The ropp_fm processing is:
ropp_fm_roprof2obs:  y%geop = ro_data%Lev2a%geop_refrac
ropp_fm_obs2roprof:  ro_data%Lev2a%geop_refrac = y%geop
                     ro_data%Lev2a%alt_refrac = geopotential2geometric()
i.e. only the geopotential heights are used in the processing, and
alt_refrac is calculated from geop for output.

The difference in alt_refrac you observe is indeed due to slightly
different methods used in the GRAS SAF NRT processing and ROPP for
converting between geometric and geopotential height scales. I have
scanned through the NRT code (function Alt_from_Geop in Lib/Earth.f90)
but it is not clear from this what the exact method used is - no
references as far as I can see. However, it is similar (but more
complicated) to that used in ROPP (maybe ROPP is a simplified version?).
The ROPP routine is in ropp_tools/geodesy/geopotential2geometric (note
ropp_tools is renamed ropp_utils in ROPP-2). More details are in GSR 02
- to appear online soon! 

I do not propose we change the ROPP-2 code ahead of release. It would be
good to understand if either ROPP or GRAS SAF NRT methods is 'best' and
make any updates for future versions of either software codes.

This ticket is a reminder to look at this issue (and potentially update either GRAS SAF NRT processing or ropp_utils for consistency).

[Huw Lewis]

Short report comparing ROPP and GRAS NRT (i.e. Invert) geodesy routines

[Huw Lewis]

Initial investigations - see attached document.

Summary:

• I think both packages are based on Somagliana's equation.

• ROPP uses a constant value for equatorial gravity of 9.7803253359 m/s2 but Invert calculates this from other parameters which give 9.7805898 m/s2. This does not look right to me - can Michael explain this?

• The different g_e values result in a constant offset between computed surface gravity values at all latitudes.

• This results in differences between converted geopotential height values of up to 1.5~m at 60 km when converted using ROPP and Invert routines.

• The simplified ROPP expressions are just as 'accurate' in terms of the resulting geop<->geom height conversions.

[Huw Lewis]

Response from Michael Gorbunov:

The approximation for the gravity field that I adopted was based on the following two sources:

1. K. Lambeck, Geophysical Geodesy, 1988 
2. US Standard Atmosphere



I also have a module that implements a more accurate gravity field based on geoid model (coefficient file
 jgm3_osu91a) that was supplied by Stephen Leroy.
 The comparison reveals small differences between the two models.

I do not have access to Lambeck (1988). Refering to W. Torge (2001) "Geodesy", published by Walter de Gruyter, give values of normal gravity at equator = 9.7803267715 m/s2. This is still in better agreement with the ROPP value. Will maintain 'ROPP' constants for ROPP-3 release, but ticket remains open until issues are clarified further.

Need improved commenting of relevant code comments with full references for all constants used.

[Huw Lewis]

At occ code v19.4.571, Michael has updated his coefficients to be in line with those used in ROPP. Thus, we expect this discrepancy should now be removed.

Ticket remains open for more thorough assessment of this difference between the two versions.

[Huw Lewis]

Testing this consitency (e.g. in output test folder) is now made more difficult by new OCC licence agreement with Michael Gorbunov, which prevents direct write-out of OCC results with ropp_io_write. Instead, we have to deal with ascii output (with height data rounded to nearest m) and convert these to ROPP format using occ2ropp tool.

The more important issue is whether the ROPP routines are sufficiently accurate. See #159.

Ticket closed as invalid.