id	summary	reporter	owner	description	type	status	priority	milestone	component	version	resolution	keywords	cc
138	Thinning refractivity bugfix	Huw Lewis	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
}}}
"	defect	closed	normal	3.0	ropp_io	1.2	fixed		dave.offiler@… kbl@… mbs@… ssy@… frr@… huw.lewis@…