DMI ROPP9.0 wishlist (BAROCLIM etc)

[Ian Culverwell]

Stig Syndergaard (DMI) writes:

Hi Ian,

Sorry for this rather lengthy email. For some time now, I wanted to 
update you on the dmi_trunk ROPP work here at DMI, but didn't want 
to spam you with every little thing. So you get everything in a big chunk! 
Hope you don't mind.


We finally got BAROCLIM into the dmi_trunk:

https://trac.romsaf.org/ropp/changeset/4527/ropp_src/branches/dev/Share/dmi_trunk_8.0

Helge did most of the work to implement BAROCLIM. On top of this I 
added a new 2-parm fit approach for the stat.opt. (the one I 
developed in the struct.uncert. project some years back). I did 
this because the current 2-parm fit approach in ROPP has some minor 
flaws (sometimes not fitting correctly because of too few points to 
fit to) that I believe would overshadow any improvement we might 
see with BAROCLIM over MSIS when using the 2-parm fit.

Both are implemented as options, so one can choose either MSIS or 
BAROCLIM (in the configuration file), in combination with old or 
new 2-parm fit (or old 1-parm fit).

For the new 2-parm fit approach there are a couple of new
subroutines:

ropp_pp_fit_model_refraction_new.f90
ropp_pp_search_model_refraction_new.f90

with '_new' following the precedence of Chris B.'s implementation of
ropp_fm_refrac_1d_new.f90 etc...

In the .cf files the new method is called 'regular' and the old one 
'convoluted'. This because the old search and fit method is in fact 
convoluted in that it is fitting within the search algorithm. The 
new approach does first the search, then the fit. I'd be happy to 
explain more details if necessary.


The work on the 2-parm fit was done here:

https://trac.romsaf.org/ropp/changeset/4526/ropp_src/branches/dev/Share/ssy_8.0

The work on BAROCLIM was done here:

https://trac.romsaf.org/ropp/changeset/4334/ropp_src/branches/dev/Share/hjs_6.1
https://trac.romsaf.org/ropp/changeset/4482/ropp_src/branches/dev/Share/hjs_6.1
https://trac.romsaf.org/ropp/changeset/4518/ropp_src/branches/dev/Share/ssy_8.0
https://trac.romsaf.org/ropp/changeset/4519/ropp_src/branches/dev/Share/ssy_8.0
https://trac.romsaf.org/ropp/changeset/4520/ropp_src/branches/dev/Share/ssy_8.0

In the process of all this I found a small bug that in some 
compilations (when I didn't compile from scratch after a merge) 
seemed to make havoc:

https://trac.romsaf.org/ropp/changeset/4528/ropp_src/branches/dev/Share/dmi_trunk_8.0

There are also a few older changes by Hans, Joe, and myself that dates back to
summer last year (but after the 'deadline' of ROPP 8.0 wishes):

Changes to Tdry calculation by Hans:

https://trac.romsaf.org/ropp/changeset/4287/ropp_src/branches/dev/Share/dmi_trunk_6.1

Changes to background covariance calculation by Joe:

https://trac.romsaf.org/ropp/changeset/4331/ropp_src/branches/dev/Share/dmi_trunk_7.1

Changes related to eum2ropp by myself:

https://trac.romsaf.org/ropp/changeset/4278/ropp_src/branches/dev/Share/dmi_trunk_6.1
https://trac.romsaf.org/ropp/changeset/4412/ropp_src/branches/dev/Share/dmi_trunk_7.1

Hope you will consider taking all this on board for ROPP 9.0. We  should discuss
at next ROPP-GG. I understand that there will also be testing and documentation
that needs to be done. Of course, we  have tested BAROCLIM and the new 2-parm
fit here at DMI with a few profiles and some stats, but I suppose you would also
want to do some final testing when implemented in ROPP 9.0. I'd be happy to help
with these things.

Seems fairy nuff.

[Ian Culverwell]

Further input from Stig, 30/7/2015:

Hi Ian,

Running BAROCLIM with ropp_pp_occ_tool here at DMI 
revealed a bug (Hans found it; I had only tested with 
ropp_pp_invert_tool).

I had forgotten that the phase model calculation under 
ropp_pp/preprocess is always based on MSIS 
refractivity, there is no equivalent for BAROCLIM (it's 
only a BA model).

I have introduced a new configuration parameter 
'msisfile' that solves the problem:

https://trac.grassaf.org/ropp/changeset/4531/ropp_src/branches/dev/Share/dmi_trunk_8.0

Sorry about this late addition.

-Stig

[Ian Culverwell]

Stig also asks, in connection with Axel's changes to eum2ropp,

BTW, I agree very much with this comment around line 
2410 in ropp_io_read_ncdf_get.f90:

'It would be better to store the reference values for 
e.g. the bufr in a dedicated dimension, not in the 
level 1a one! This dimension would always have the 
value 1 when used, but allows to separate data at the 
reference point from the general level 1a data, which 
covers the whole occultation.'.

Are you considering that for ROPP 9.0?

We can certainly consider it, but I'm not sure how much code engineering would be needed.

[Stig Syndergaard]

File with open loop data for testing eum2ropp

[Ian Culverwell]

Hans' Tdry calculation included at r4781.

[Ian Culverwell]

The code to extend the -l option to eum2ropp (basically Stig's r4278) has been merged into the ROPP90_prototype branch at r4825.

[Ian Culverwell]

The appropriate default for the -l option of eum2ropp remains to be decided. Currently it's none, which has the same effect as not using -l, and which is different to the pre-ROPP9.0 default, which was equivalent to the new -l cl+rs. To be discussed with Stig and Axel (the only likely users of this tool).

[Ian Culverwell]

lev1a_plot_p1.gif

[Ian Culverwell]

lev1a_plot_p2.gif

[Ian Culverwell]

The phase_L1 produced by eum2ropp -l cl, eum2ropp -l cl+ol and eum2ropp -l cl+rs on Stig's example file (which has cl, ol and rs data in it) is

Looks sensible: no nasty jumps across the 'type' boundaries.

Similarly for phase_L2:

[Ian Culverwell]

The BAROCLIM changes described above (r4527 + r4528 + r4531) have been merged into the ROPP9.0 prototype branch at r4830.

[Ian Culverwell]

Testing of the BAROCLIM changes revelaed one anomaly. When using -m MSIS (ie, a local search), the new code reads the bending angles from the MSIS_coeffs.nc file, but the old code read refractivities from this file, interpolated them to high vertical resolution, and calculated the bending angles from them using the Abel transform. These two methods give slightly different results. It is not yet clear whether we live with this change to the default behaviour, or if the means are provided to allow the user to revert to the old behaviour, if desired.

[Ian Culverwell]

Documentation revisions have been committed at r4833.

[Ian Culverwell]

Further minor documentation changes have been committed at r4834, r4835 and r4836.

[Ian Culverwell]

ROPP9.0 PP UG

[Ian Culverwell]

ROPP9.0-ROPP8.0 PP UG

[Ian Culverwell]

ROPP8.1 PP UG

[Ian Culverwell]

latexdiff.sh

[Ian Culverwell]

The (provisional) ROPP9.0 PP UG, the ROPP8.1 PP UG, and the differences between them (according to latexdiff), ROPP9.0-ROPP8.1 PP UG, are attached.

Note that latexdiff doesn't pick up differences in verbatim latex. Also note that Fig 4.2 has been updated ('MSIS' --> 'climatological').

[Ian Culverwell]

For the record, here are some differences between the ROPP9.0 implementation

ropp_pp_invert_tool  IT-PP-02_prof1.nc  -o IT-PP-02_prof1_test2.nc  -mfile BAROCLIM_coeff.nc  -c grassaf_invert_test2.cf  -m BARO

(ie the BAROCLIM options, with sf_method = regular), and the earlier, ROPP8.1,

ropp_pp_invert_tool  IT-PP-02_prof1.nc  -o IT-PP-02_prof1_cntl.nc  -mfile MSIS_coeff.nc  -c grassaf_invert_cntl.cf  -m MSIS

(ie the MSIS options, with sf_method = convoluted).

[Ian Culverwell]

ropp_pp_invert_tool_BAROCLIM-MSIS.png

[Ian Culverwell]

IT-PP-02_prof1.nc

[Ian Culverwell]

IT-PP-02_prof1_cntl.nc

[Ian Culverwell]

IT-PP-02_prof1_test2.nc

[Ian Culverwell]

grassaf_invert_cntl.cf

[Ian Culverwell]

grassaf_invert_test2.cf

[Ian Culverwell]

And here are the corresponding differences between the BAROCLIM options of ropp_pp_occ_tool,

ropp_pp_occ_tool  IT-PP-04_prof1.nc  -o IT-PP-04_prof1_test2.nc  -mfile BAROCLIM_coeff.nc  -c cosmic_pp_test2.cf  -ellipsoid  -m GBARO

and the MSIS options,

ropp_pp_occ_tool  IT-PP-04_prof1.nc  -o IT-PP-04_prof1_cntl.nc  -mfile MSIS_coeff.nc  -c cosmic_pp_cntl.cf  -ellipsoid  -m GMSIS

[Ian Culverwell]

IT-PP-04_prof1.nc

[Ian Culverwell]

IT-PP-04_prof1_cntl.nc

[Ian Culverwell]

IT-PP-04_prof1_test2.nc

[Ian Culverwell]

cosmic_pp_cntl.cf

[Ian Culverwell]

cosmic_pp_test2.cf

[Ian Culverwell]

ropp_pp_occ_tool_BAROCLIM-MSIS.png

[Ian Culverwell]

So the ropp_pp_occ_tool changes are around 10 times larger than those from ropp_pp_invert_tool.

[Ian Culverwell]

Note that there's also a difference between ROPP8.1 and ROPP9.0 when running ropp_pp_invert_tool (and, in fact, ropp_pp_occ_tool) with the default options (but with -m MSIS):

This is a result of the different methods of generating bending angle climatologies at ROPP8.1 (calculated by Abel inversion of the climatological refractivities) and ROPP9.0 (calculated directly from the coefficients in the file). As can be seen, the differences are around 10 times smaller than the BAROCLIM-MSIS differences shown above.

[Ian Culverwell]

ropp_pp_invert_tool_ROPP9.0-ROPP8.1.png

[Ian Culverwell]

(Note that there is no effect, on either ropp_pp_invert_tool or ropp_pp_occ_tool, if one uses -m GMSIS or -m GMBARO, because global searches read the bending angle directly from the climatological file, and this hasn't changed between ROPP8.1 and ROPP9.0.)

[Ian Culverwell]

Updates to documentation lodged at r4864, and to ropp_pp_occ_tool help information and man page at r4865.

[Ian Culverwell]

Raised possibility of removing the sf_method = convoluted option from ROPP, since the -m MSIS option has broken backwards compatibility anyway (see above). Since I've just gone to the trouble of documenting it, I think I'll probably leave it in - at ROPP9.0, at least. To be decided.

[Ian Culverwell]

Some users may require the old method, so decide, for the moment to leave both fitting methods in place. But, after some discussion it was agreed to change the names of the two methods, so that 'regular' is now called 'serial', and 'convoluted' is now called 'parallel'. The latter remains the default.

Code changes check out OK, so they have been committed at r4931. The ROPP PP UG changes have been made at r4932.

[Ian Culverwell]

BAROCLIM3 licence added to distro COPYRIGHT files at r4957.

BAROCLIM3 licence added to ROPP documentation at r4958.

[Ian Culverwell]

It turns out that the BAROCLIM_coeffs.nc I was using was version 2, not version 3. The later version has been committed to the repository at r4959.

[Ian Culverwell]

Results are pretty similar to BAROCLIM (v2). For ropp_pp_io_invert(BAROCLIM - MSIS) we now have:

Same smallish differences, except for Tdry, but that's not my problem.

Similarly for ropp_pp_io_occ(BAROCLIM - MSIS) we now have

Ditto.

[Ian Culverwell]

ropp_pp_occ_tool_BAROCLIM-MSIS2.png

[Ian Culverwell]

ropp_pp_invert_tool_BAROCLIM-MSIS2.png

[Ian Culverwell]

(I meant ropp_pp_occ_tool, of course, not ropp_pp_io_occ. Sim. for invert.)

[Ian Culverwell]

It turns out the r4521 was accidentally included in the BAROCLIM3 build, so remove that at r4985.

[Ian Culverwell]

These requests have now all been included in the ROPP90 prototype, so this ticket is being closed as 'fixed'.

[Ian Culverwell]

Hans has also made some changes atmPrf/atmPhs netCDF readers, to generate more representative ROPP PCDs. These changes have been made at r4930. (Note that these changes also include a correction to the value of g:

g_wmo = 9.80665_wp

rather than

g_wmo = 9.80665

This could produce ~1e-07 fractional differences in the results.

[Ian Culverwell]

Reclosing ticket as fixed.

[Ian Culverwell]

Reopening to log two minor bug fix requests from Stig:

I should notify you also of a few bugs we found in the 
OCC code when we did the WO data validation. We are 
still using Michael's OCC code in NRT, but the same 
bugs are in ROPP:

https://trac.romsaf.org/ropp/changeset/4775/ropp_src/branches/dev/Share/dmi_trunk_8.1/ropp_pp/icorr/ropp_pp_ionospheric_correction.f90

https://trac.romsaf.org/ropp/changeset/5065/ropp_src/branches/dev/Share/dmi_trunk_8.1/ropp_pp/tools/ropp_pp_occ_tool.f90

This latter bug has only been fixed in the occ_tool in 
our dmi_trunk, since we stopped maintaining the 
invert_tool after combining invert and occ tools (using 
only the occ_tool now). It should be fairly easy to 
find the corresponding bug in ropp_pp_invert_tool.f90.

We'll try and slip these into ROPP9.0.

[Ian Culverwell]

Note that r4775 would only come into force if config%dpi were >= 50, and it equals 100 in all our standard config files, so the change should never be invoked in our tests. But if we arbitrarily set it to 25 (m) to test things, then ropp_pp_occ_tool runs OK and the resulting differences in bending angle

Image(bangle_dpi100_and_dpi25.png)

and refractivity

Image(refrac_dpi100_and_dpi25.png)

look sensible.

[Ian Culverwell]

I'll try all that again! (Really a shame we can't edit tickets post-cock-ups.)

Note that r4775 would only come into force if config%dpi were <= 50, and it equals 100 in all our standard config files, so the change should never be invoked in our tests. But if we arbitrarily set it to 25 (m) to test things, then ropp_pp_occ_tool runs OK and the resulting differences in bending angle

and refractivity

look sensible.

[Ian Culverwell]

bangle_dpi100_and_dpi25.png

[Ian Culverwell]

refrac_dpi100_and_dpi25.png

[Ian Culverwell]

refrac_dpi100_cf_dpi25.png

[Ian Culverwell]

On closer inspection, there is a systematic difference in the refracs above about 30 km. (The bending angles are OK.) I don't know what this means.

.

[Ian Culverwell]

Similarly little impact on invert bending angles (not plotted), and a much smaller one on the refracs:

I don't know what this means, either.

[Ian Culverwell]

refrac_dpi100_cf_dpi25_invert.png

[Ian Culverwell]

As Stig rightly suggests, this is because I had to use GMSIS to avoid the ROPP8.1-ROPP9.0 difference I discussed above. This results in different climatologies appearing in the different searches. When I do it properly, the refractivity differences in occ

and invert

largely disappear. (Note much expended difference scales.)

[Ian Culverwell]

refrac_dpi100_cf_dpi25_occ.png

[Ian Culverwell]

refrac_dpi100_cf_dpi25_invert2.png

[Ian Culverwell]

Change committed at r5104.

[Ian Culverwell]

Now look at r5065. The change from

            WHERE ( out_ba%impact_L1 > config%Pmax - 5000.0_wp )

to

            WHERE ( out_ba%impact_L1 > out_ba%impact_L1(imax) )

in ropp_pp_occ_tool.f90, means, for IT-PP-04_prof1.nc,

 config%Pmax - 5000.0_wp  =   6516374.58453080     
 out_ba%impact_L1(imax)   =   6516387.27594711     
 MAXVAL(out_ba%impact_L1) =   6544074.50490837     
 COUNT(out_ba%impact_L1 > config%Pmax - 5000.0_wp) = 278
 COUNT(out_ba%impact_L1 > out_ba%impact_L1(imax))  = 277

The change therefore makes ~1e-4 difference to the last point in bangle_L{1,2}, and ~1e-15 diff to the last few points of bangle_opt:

There's no change in the refracs.

The bugfix has a similar or smaller impact on the other 19 profiles in IT-PP-04.nc. So it seems a pretty safe change.

However, however, however, for ropp_pp_invert_tool there are much larger, more widespread differences between bangle_opt

and refrac

This is the case even though the differences appear to be restricted to the last point in the bangle_opt profile, as for ropp_pp_occ_tool:

 config%Pmax - 5000.0_wp  =    6432757.62800000     
 out_ba%impact_L1(imax)   =    6432795.71606681     
 MAXVAL(out_ba%impact_L1) =    6503967.86234704     
 COUNT(out_ba%impact_L1 > config%Pmax - 5000.0_wp) = 713
 COUNT(out_ba%impact_L1 > out_ba%impact_L1(imax)) =  712

Oh blimus.

[Ian Culverwell]

r5056_bangle_opt_occ.png

[Ian Culverwell]

r5056_bangle_opt_invert.png

[Ian Culverwell]

r5056_refrac_invert.png

[Ian Culverwell]

I think it may be connected with the fact that, for IT-PP-02_prof1.nc,

imin, imax =            1         720
COUNT(out_ba%impact_L1 > config%Pmax - 5000.0_wp) =  713
COUNT(out_ba%impact_L1 > out_ba%impact_L1(imax))  =  712

In other words, almost all the L1s (and L2s) are being set to the model values anyway. This doesn't sound right. (For occ, the figures were

imin, imax =            1        1190
COUNT(out_ba%impact_L1 > config%Pmax - 5000.0_wp) = 278
COUNT(out_ba%impact_L1 > out_ba%impact_L1(imax))  = 277

which seems more sensible.)

So I tried using the output of OCC (1190 dim_lev1b points) as the input to INVERT. I find

imin, imax =            1        1140
config%Pmax - 5000.0_wp =    6511387.27594711     
out_ba%impact_L1(imax) =    6511389.58118876     
MAXVAL(out_ba%impact_L1) =    6544074.50490837     
COUNT(out_ba%impact_L1 > config%Pmax - 5000.0_wp) = 328
COUNT(out_ba%impact_L1 > out_ba%impact_L1(imax))  =  327

which looks better. And indeed, there are much smaller errors in bangle_opt

and refrac

are now much smaller. But it's all a bit strange.

[Ian Culverwell]

r5056_bangle_opt_invert2.png

[Ian Culverwell]

r5056_refrac_invert2.png

[Ian Culverwell]

I was wrong with my accounting above: there are nbi=1432 points in the profile, and imax is the point closest to config%Pmax-5000.0_wp, so what we're actually saying is that about half of the points are set equal to the model background, not almost (or more than!) all of them.

Stig said

Regarding the bigger differences for invert_tool than 
for occ_tool when fixing the bug: It could be the 
number of points that makes the difference. The 
occ_tool results in a denser grid (for whatever reason) 
with 1190 points from imin to imax instead of 720. That 
could maybe explain why the effect is smaller for the 
occ_tool and when running invert_tool with the denser 
grid. Not sure, but it seems plausible.

I replied

I agree that the different density of data points might 
explain the different behaviours of occ and invert when 
the bug is fixed or not.  It puzzles me that the 
changes in refracs are _so_ much bigger for invert 
(~1e-3) than occ (~1e-10), but there we go.

Stig replied

hmm... yes, I agree, something there is not understood. 
Also that you see such big difference further down the 
profile (from your plots in the ticket). Perhaps if an 
integer number of points in some filtering is changed 
by 1, and if that happens more often in invert_tool 
because of the lower density of points.... I don't know.

These large changes mean that we would probably need to make significant changes to the reference files in the Test Folder, and it's a bit late to do that for ROPP9.0.

The bugfix is certainly a sensible thing to do, however. The attached pfixcases.ps file, pinned to this ticket with Stig's permission, shows the impact of the bug in 7 cases where the bug had such a large influence that the refractivity became negative at the top sample (Stig found only a few such cases over two days that he looked at - 4 on 20160809 and 3 on 20160806). Each case consists of three plots:

'....XXXX...' in legends came from DMI's production (or WO test)

'....ROPP...' in legends came from ROPP before the bugfix

'....PFIX...' in legends came from ROPP after the bugfix.

Red: WO

Green: GO NRT for reference

The bugfix is certainly beneficial, although it doesn't completely cure the negative refracs (eg profs 3, 4 and especially 5). But even in these cases it's much better than the buggy ROPP. So we should definitely try to include it.

Conclusion: Defer this necessary bugfix to ROPP10.0.

[Ian Culverwell]

Moving ticket to ROPP10.0.

[Ian Culverwell]

pfixcases.ps

[Ian Culverwell]

Changes to the devt branch have been made at r5595 and r5596. (The latter contains some polishing of the testing routines too.)

[Ian Culverwell]

I don't think this caused any issues with the test folder, or if it did they've been updated, so close the ticket.