Thinning refractivity bugfix

[Huw Lewis]

Email from Stig

We may have found a problem in ropp_io_thin.f90 regarding the thinning 
of refractivity. Both Frans, Kent and Martin has contributed to finding 
this problem. I spend most of the day looking in part of the code and 
here is my assessment:

Under point 3 in the code there is the following:

! Actually we thin on Impact Parameter levels (ThinLev left from BA 
thinning)
! as this is independent of Level 1b profiles; converting ThinLev to
! geometric height here would need Impact Parameters which may not exist in
! the ROdata structure, or could be on different levels to Refractivity 
levels.
! After thinning, Impact Parameter levels (still in ThinLev) can then be
! converted to geometric height using the thinned refractivites.
!
! Convert full refractivity levels in geometric ht (wrt geoid) to IP using:
!  1) h = r-Rc
!  2) a = n.r
!  3) n = 1+N.10^-6
!  => a = (1+N.10^-6).(h+Rc)
!
! If N is invalid (which may be the case at high and low fixed levels),
! use a log interpolation between height of last good value (initialised
! with N=10^-6 at highest level) and a reference value (N=300) at the 
surface
! (h=0). This is just to give sensible impact parameter values for all of
! the fixed levels which can be inverted back to geometric heights later;
! they do not need to be accurate, since by definition, the refractvity is
! missing so these levels won't be useful anyway.


And the coding right below is according to these equation. But is there 
not a mistake there? h = r-Rc is wrt ellipsoid, not geoid. In the code 
there is a line:

        Lev(i) = ( 1.0_wp + 1e-6_wp * N ) &
                * ( ROdata%Lev2a%Alt_refrac(i) + ROdata%GEOref%RoC )

But since Alt_refrac is wrt to the geoid (I assume it is the variable 
coming from the NetCDF file which is now wrt geoid), Lev is then not the 
impact parameter but will be offset by the undulation (approximately). 
The thinning that follows will then be on wrong levels. Later in the 
code (after the thinning) the correct equations are used (including the 
undulation):

         ThinLev(i) = ( ThinLev(i) / ( 1.0_wp + 1e-6_wp * N ) ) &
                    - ROdata%GEOref%RoC &
                    - ROdata%GEORef%Undulation

I think the net result is that the thinned profile is offset by the 
undulation (approximately) as compared to the hi-res profile.

Should the undulation not be included in the first part as well?


Unrelated to this (and less critical, but still something that we should 
fix in some next version of ROPP), both of the do loops where the above 
lines are found assumes that N = 300 when H = 0. This is used for 
interpolation/extrapolation of refractivity to find levels where there 
is no data. I understand that these N values are used only for getting 
fairly realistic values of all 247 levels, even when data are missing, 
okay! But the H here is impact height, not altitude; as you know there 
is a difference of a few km near the surface. So the N near the surface 
will be largely underestimated. It will actually be more like 225 if the 
last good N was at a high altitude (assuming 2 km difference between 
impact height and altitude). This results in about 560 km for the lowest 
level in the BUFR, not near zero as would be more reasonable. Frans has 
observed similar values for the lowest level in the BUFR files when data 
are missing. Perhaps N=400 at H=0 in the code would be more in line with 
reality (for the purpose of interpolation to values of H ~2km and above) 
when H is impact height.

Regards,
-Stig

[Huw Lewis]

All,

Undulation correction in refractivity thinning code

Having looked through the ropp_io_thin code, I agree with Stig's thoughts that the alt_refrac data need to have the undulation subtracted before calculating impact parameters. This has been modified in the code, so the calculations of impact parameter from alt_refrac and later alt_refrac from thinned impact parameters are consistent. See revision ​https://trac.grassaf.org/ropp/changeset/1734/ropp_src/branches/dev/Share/r1661_ropp_pp/ropp_io/thin/ropp_io_thin.f90.

It is not clear why only one of these calculations was modified when the thinner routines were first written. This revision [1105] shows when undulation was applied to the calculation of alt_refrac after thinning "Clarified that alt_refrac parameter should be wrt geoid, not ellipsoid" ​https://trac.grassaf.org/ropp/changeset/1105/ropp_src/trunk/ropp_io/ropp/ropp_io_thin.f90 Both calculations were consistent before this change. It APPEARS that it was realised that the need to compute alt_refrac wrt geoid resulted in the code change, but this was only applied in one place. Perhaps at this time it was not clear how the incoming alt_refrac data were defined? Ticket remains open until further testing with new solution found to produce more consistent results.

Use of N=300 at surface for interpolating where missing data

This code has also been reviewed. As stated, the code is considering N=300 at an impact height of 0 (rather than at the surface where geometric height=0). While this is going to lead to numerical errors in height, the height values at these locations is not thought to be important for data quality, since no refractivity data are available at these levels.

Code review has identified that we are actually having to use impact parameter values from Level1b data here in any case, which may not correspond with the levels in the Level2a data (or even be ordered in the same way). The assumption that there are at least as many Level1b points as Level2a required for this code to work is considered reasonable. Ticket remains open as a reminder to revisit the thinner code and review in more detail.

[(none)]

Milestone 3.0 deleted

[Huw Lewis]

Comparison of the MT-IO-03 thinner test between v1.2 and v2.0 show the impact of the changes to the ROPP thinner between these versions. See files attached to #134 for plots.

The pre-sorting of impact parameters before thinning reduces the standard deviation of thinned profiles for bending angle with most thinning methods. The largest difference between v1.2 and v2.0 results from correcting the undulation bugfix in thinning refractivity (note the scale on the axes reduce from 0.6 to 0.1 between plots for v1.2 and v2.0). There is also a small change in the distribution of standard deviations for thinned temperature profiles, but note larger biases at low levels for v1.2 are largely removed in v2.0 comparisons.

Ticket closed as fixed - will open a new ticket to review thinner code and produce GSR.