Incorporation of ropp_pp_invert_tool in ropp_pp_occ_tool [5 d]

[Ian Culverwell]

https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_pp/common/Makefile.am https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_pp/common/ropp_pp.f90 https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_pp/common/ropp_pp_diag2roprof.f90 https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_pp/common/ropp_pp_interpolate_latlonaz.f90 https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_pp/common/ropp_pp_read_config.f90 https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_pp/common/ropp_pp_types.f90 https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_pp/icorr/ropp_pp_fit_model_refraction.f90 https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_pp/icorr/ropp_pp_fit_model_refraction_new.f90 https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_pp/icorr/ropp_pp_ionospheric_correction.f90 https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_pp/tools/ropp_pp_invert_tool.f90 https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_pp/tools/ropp_pp_occ_tool.f90

Original changesets:

Incorporation of ropp_pp_invert_tool in ropp_pp_occ_tool: https://trac.romsaf.org/ropp/changeset/4638/ropp_src/branches/dev/Share/dmi_trunk_8.0/ropp_pp/tools/ropp_pp_occ_tool.f90

[Ian Culverwell]

See r6788 (mainly), r6791 and r6792.

[Ian Culverwell]

Introduce a new test script, test_pp_occ_inv.sh, which tests the -inv option to ropp_pp_occ_tool, at r6848.

[Ian Culverwell]

The outputs of ropp_pp_invert_tool and ropp_pp_occ_tool -inv are identical apart from the following.

1) occ -inv includes the extra variables

	float IC_badness_score(dim_unlim) ;
		IC_badness_score:long_name = "Ionospheric correction badness score" ;
		IC_badness_score:units = "1" ;
		IC_badness_score:valid_range = 0., 1000. ;
	float impact_sobgr(dim_unlim, dim_lev1b) ;
		impact_sobgr:long_name = "SO background impact parameters" ;
		impact_sobgr:units = "metres" ;
		impact_sobgr:valid_range = 6200000., 6600000. ;
	float bangle_sobgr(dim_unlim, dim_lev1b) ;
		bangle_sobgr:long_name = "SO background bending angle" ;
		bangle_sobgr:units = "radians" ;
		bangle_sobgr:valid_range = -0.001, 0.1 ;

This is expected, I think.

2) The global attributes of invert include

		:bangle_method = "UNKNOWN" ;
		:refrac_method = "UNKNOWN" ;

whereas the occ -inv equivalents are

		:bangle_method = "UNKNOWN (StatOpt)" ;
		:refrac_method = "Abel transform (LIN)" ;

Again, probably expected or acceptable.

3) Significant differences between lon_tp and lat_tp: invert gives

    lat_tp = 
    -9.9999e+07, -9.9999e+07, -9.9999e+07, -9.9999e+07, -9.9999e+07, ...

    lon_tp = 
    -9.9999e+07, -9.9999e+07, -9.9999e+07, -9.9999e+07, -9.9999e+07, ...

whereas occ -inv gives

    lat_tp = 
    -1.42127e-09, -1.42127e-09, -1.42127e-09, -1.42127e-09, -1.42127e-09, ...
    lon_tp = 
    -4.9032e-09, -4.9032e-09, -4.9032e-09, -4.9032e-09, -4.9032e-09, ...

(Note that

    lat_tp = 
    -1e+08, -1e+08, -1e+08, -1e+08, -1e+08, ...
    lon_tp = 
    -1e+08, -1e+08, -1e+08, -1e+08, -1e+08, ...

in the input file.)

I think this needs fixing, or addressing at least.

[Ian Culverwell]

Input file ropp_pp_test.nc attached.

invert output file ropp_pp_test_1m.nc attached.

occ -inv output file ropp_pp_test_1n.nc attached.

(Default options in both cases: ropp_pp_invert_tool ropp_pp_test.nc -o ropp_pp_test_1m.nc -d and ropp_pp_occ_tool ropp_pp_test.nc -inv -o ropp_pp_test_1n.nc -d)

[Ian Culverwell]

ropp_pp_test.nc

[Ian Culverwell]

ropp_pp_test_1m.nc

[Ian Culverwell]

ropp_pp_test_1n.nc

[Ian Culverwell]

(Trivially, invert says

INFO:  Retrieving bending angle profile by STATISTICAL OPTIMISATION 

whereas occ -inv says

INFO:  Correcting bending angle profile by STATISTICAL OPTIMI
SATION 

There's also a bit more info in the occ -inv output - see the attached log files t_pp_invert_1.log and t_pp_occ_inv_1.log.)

[Ian Culverwell]

t_pp_invert_1.log

[Ian Culverwell]

t_pp_occ_inv_1.log

[Ian Culverwell]

All seems to pan out OK in practice, so closing ticket. (Don't know why this wasn't done earlier.)

[Ian Culverwell]

Reopening ticket as Stig needs to look at the problem with lon_tp and lat_tp that was identified in #comment:4.

[Ian Culverwell]

Stig happy to leave bangle_sobgr and impact_sobgr in ROPP-11 (for now). They were removed from DMI's ropp_pp_occ_tool.f90 at r4925, where the outputting of SO bangle was transferred to ropp_pp_diag2prof.f90. (We have the same development in ROPP-11, but it isn't activated.)

[Ian Culverwell]

The root of the crazy lon_tp and lat_tp is the new interpolation routine ropp_pp_interpol_latlonaz.f90. These fields are missing (i.e. equal to ropp_MDFV = -99999000.000) in the input file. In earlier versions of ropp_pp_occ, these were individually interpolated vertically from the old to the new impact parameter grids using ropp_pp_interpol. But since they were constant in the vertical, this interpolation had no effect, and the resulting lon_tp and lat_tp were still ropp_MDFV.

But ropp_pp_interpol_latlonaz.f90 combines the fields before interpolating, to convert to cartesian coordinates, and the numerical imprecision of doing this with angles ~ -10⁸ deg causes the funny values.

Specifically, the (undocumented and apparently unreviewed) ropp_pp_interpol_latlonaz.f90 code says

  DO i=1,nx
   rlat = Pi/2.0_wp - lattp(i)*deg2rad   ! 
   rlon = lontp(i)*deg2rad
   uvec(1,i) = sin(rlat)*cos(rlon)
   uvec(2,i) = sin(rlat)*sin(rlon)
   uvec(3,i) = cos(rlat)

   IF(use_azm) THEN
    razm = azmtp(i)*deg2rad
    vvec(1,i) = -1.0_wp*cos(rlat)*cos(rlon)*cos(razm) - sin(rlon)*sin(razm)
    vvec(2,i) = -1.0_wp*cos(rlat)*sin(rlon)*cos(razm) + cos(rlon)*sin(razm)
    vvec(3,i) = sin(rlat)*cos(razm)
   ENDIF

So if lon = lat = azi = ropp_MDFV = -99999000 = -27775 * 360 (!) then we should have

u = (1, 0, 0)

and

v = (0, 0, 1).

In practice, because of rounding error, we find

u = (1.0000000000000000, -8.5576981117344849E-011, -2.4805976052282926E-011)

and

v = (2.4805976044959505E-011, -8.5576981119467666E-011, 1.0000000000000000)

These are the same at every level, and are therefore returned unaltered by the ropp_pp_interpol interpolation routine.

But the inversion of these u and v to lon, lat and azi returns

lon_int(1:5) = -4.9031998414945999E-009 -4.9031998414945999E-009 -4.9031998414945999E-009 -4.9031998414945999E-009 -4.9031998414945999E-009

lat_int(1:5) = -1.4212702126314980E-009 -1.4212702126314980E-009 -1.4212702126314980E-009 -1.4212702126314980E-009 -1.4212702126314980E-009

azm_int(1:5) = 359.99999999509680 359.99999999509680 359.99999999509680 359.99999999509680 359.99999999509680

rather than ropp_MDFV as we would want.

(Note that even if the code did the calculation and inversion perfectly, it would still deliver lon = lat = azi = 0 (mod 2 pi), because the input has lon = lat = azi = 0 (mod 2 pi).)

The solution is clearly to mask out the missing data - something that was not required before, for the reasons discussed above.

[Ian Culverwell]

Note that the azimuths of 359.99999999509680 are reset to ropp_MDFV in the output file, because the valid range of azimuth_tp in the input file, ropp_pp/data/ropp_pp_test.nc, is (0, 359.9):

	float azimuth_tp(dim_unlim, dim_lev1b) ;
		azimuth_tp:long_name = "GNSS->LEO line of sight angles (from True North) for tangent points" ;
		azimuth_tp:units = "degrees" ;
		azimuth_tp:valid_range = 0., 359.9 ;

I don't know why that is - the usual ROPP valid range of azimuths is (9, 360), e.g. see ropp_io/ropp/ropp_io_types.f90:

    REAL(wp), DIMENSION(2)          :: azimuth_tp      = (/   0.0_wp,  360.0_wp/)

Just an old file, perhaps. Anyway, it's enough to cause azimuths of 359.99999999509680 to be nullified if they are output with range-checking.

[Ian Culverwell]

Some valid data checks have been introduced to ropp_pp_interpolate_latlonaz.f90 at r6898. It checks out OK:

input: (lon, lat, azi) = (mdi, mdi, mdi) ==> output (lon, lat, azi) = (mdi, mdi, mdi) (all missing)
input: (lon, lat, azi) = (mdi, mdi,  45) ==> output (lon, lat, azi) = (mdi, mdi, mdi) (cannot process azi if lon or lat missing)
input: (lon, lat, azi) = ( 30, -60, mdi) ==> output (lon, lat, azi) = ( 30, -60, mdi) (cannot process azi if azi missing)
input: (lon, lat, azi) = ( 390, -60, 45) ==> output (lon, lat, azi) = (mdi, mdi, mdi) (input lon out of range)
input: (lon, lat, azi) = ( 30, -120, 45) ==> output (lon, lat, azi) = (mdi, mdi, mdi) (input lat out of range)
input: (lon, lat, azi) = ( 30, -60, -45) ==> output (lon, lat, azi) = (30, -60, mdi)  (input azi out of range)
input: (lon, lat, azi) = ( 30, -60,  45) ==> output (lon, lat, azi) = (30, -60,  45)  (all valid)

I'm happy that the problem I identified in comment:4 has been solved, but I'll leave the ticket open as more problems may emerge.

[Ian Culverwell]

This has been included in the ROPP11.1 prototype at r7147. Closing ticket.