Include Sentinel-6 reading code

[Ian Culverwell]

At the ROPP-11.0 DRR we agreed to include mods to allow ROPP to read Sentinel-6 data in a ROPP-11.1 development branch (not a release). This ticket will record work done for this.

[Ian Culverwell]

Following a discussion with Stig on 08/02/2022, it seems that the key changes are included in r6738 and r6895. (Some of the latter changes supersede some of the former ones.) Look into the issue of reading string arrays from netCDF files in ROPP. Check that unit_convert can work with frequencies. Think about pmode, but probably keep the switches for now. SX6 only has OL data, but treat it as CL. There's an example config file in r6831.

I should also investigate the strange behaviour of missing data in ropp_pp_occ, as highlighted in r7035.

Also try to include ropp_pp_kappa_residual. Will need documenting. Put switch in config file?

[Ian Culverwell]

1st step: download a pseudo SX6 dataset (COSMIC data with SX6 metadata) from #695. S6A_RO_1B_BND_____20200111T011620_20200111T011809_20201202T112537_R03__NN__________EUM__VAL_NT_R01.nc (attached) should do the trick.

[Ian Culverwell]

/net/home/h04/idculv/Example_RO_data/EUM/Yaros/S6A_RO_1B_BND_20200111T011620_20200111T011809_20201202T112537_R03NNEUMVAL_NT_R01.nc

[Ian Culverwell]

Next, make a branch to work on: https://trac.romsaf.org/ropp/browser/ropp_src/branches/dev/Share/ROPP111_SX6_read.

[Ian Culverwell]

First big job is to extend the ROprof structure to include signal info. (Since SX-6 (and later) satellites will be able to receive non-GPS signals, we need to read frequencies from the input file rather than assume default L1 and L2 values.) Stig has done this through a SIGtype structure, which is included in ROprof. He hasn't included frequency ranges, so range-checking won't work. I need to do that.

I also checked whether unit_convert could convert frequencies. It couldn't, so I have included that at r7092.

[Ian Culverwell]

Had a first crack at this at r7093. So far, this is just about extending ROprof to include the signal structures, updating the range-checking and ropp2ropp assignment and equivalencing routines, extending test2ropp to generate signal data (and bug-fixing it), and checking that we can read/write signal data from/to ROPP format netCDF files. It seems to work OK: test2ropp says

INFO (from test2ropp):  Generating VALID data...
 VALID: ROdata%Signal1%Range%Freq(mn), ROdata%Signal1%Range%Freq(mx), randf(41,1), ROdata%Signal1%Freq =   0.50000000000000000        5.0000000000000000       0.47722280470332246        2.6475026211649508     
 VALID: ROdata%Signal2%Range%Freq(mn), ROdata%Signal2%Range%Freq(mx), randf(42,1), ROdata%Signal2%Freq =   0.50000000000000000        5.0000000000000000       0.49748299765026782        2.7386734894262053     

which appears in the output file as

idculv@vld189:> ncdump -h ropp_test_2v.nc |grep freq
	double freq_L1(dim_unlim) ;
		freq_L1:long_name = "L1 carrier frequency" ;
		freq_L1:units = "Hz" ;
		freq_L1:valid_range = 500000000., 5000000000. ;
	double freq_L2(dim_unlim) ;
		freq_L2:long_name = "L2 carrier frequency" ;
		freq_L2:units = "Hz" ;
		freq_L2:valid_range = 500000000., 5000000000. ;

idculv@vld189:> ncks -H -Q -vfreq_L1,freq_L2 ropp_test_2v.nc
netcdf ropp_test_2v {
  dimensions:
    dim_unlim = UNLIMITED ; // (1 currently)

  variables:
    double freq_L1(dim_unlim) ;

    double freq_L2(dim_unlim) ;

  data:
    freq_L1 = 2647502621.16495 ;

    freq_L2 = 2738673489.42621 ;

The range-checking seems to be working too: test2ropp says

INFO (from test2ropp):  Generating INVALID data...
 INVALID: ROdata%Signal1%Range%Freq(mn), ROdata%Signal1%Range%Freq(mx), v, ROdata%Signal1%Freq =   0.50000000000000000        5.0000000000000000       -2084495180.5196366       -1.5844951815196366     
 INVALID: ROdata%Signal2%Range%Freq(mn), ROdata%Signal2%Range%Freq(mx), v, ROdata%Signal2%Freq =   0.50000000000000000        5.0000000000000000        439255738.99384010        5.4392557399938406     

which appears in the output file as

idculv@vld189:> ncdump -h ropp_test_1i.nc |grep freq
	double freq_L1(dim_unlim) ;
		freq_L1:long_name = "L1 carrier frequency" ;
		freq_L1:units = "Hz" ;
		freq_L1:valid_range = 500000000., 5000000000. ;
	double freq_L2(dim_unlim) ;
		freq_L2:long_name = "L2 carrier frequency" ;
		freq_L2:units = "Hz" ;
		freq_L2:valid_range = 500000000., 5000000000. ;

idculv@vld189:> ncks -H -Q -vfreq_L1,freq_L2 ropp_test_1i.nc
netcdf ropp_test_1i {
  dimensions:
    dim_unlim = UNLIMITED ; // (1 currently)

  variables:
    double freq_L1(dim_unlim) ;

    double freq_L2(dim_unlim) ;

  data:
    freq_L1 = -99999000 ;

    freq_L2 = -99999000 ;

[Ian Culverwell]

First changes to ropp_io_read_ncdf_get.f90 included in r6738 have been included in r7094. This still cannot read my example SX-6 file because

---------------------------------------------------------------------
                  EUMETSAT to ROPP netCDF Converter
---------------------------------------------------------------------

INFO (from eum2ropp):  Reading file ../data/sx6_test.n4
INFO (from ropp_io_read_eumdata):  Format Version 14.0
 pmode = NRT                           
ERROR: Variable not found: /quality/cl_data_available

This might be remedied by the changes in r6895, which includes lines like

This should be fixed by r6895, which has lines like
    IF ( ncdf_isvar('/quality/cl_data_available') ) THEN
      CALL ncdf_getvar('/quality/cl_data_available', readbyte1)
   ELSE IF (data%leo_id(1:3) == 'SE6') THEN
...

Note that I have replaced

    frequencies = (/1575420000_wp, 1227600000_wp/)

by

    frequencies = (/1.57542E9_wp, 1.2276E9_wp/)

which is closer to the ROPP standard, and makes no difference (I checked with the attached test_spdp.f90).

There was also an oddity with ol_data_used vs ol_data_available. I had to replace

      IF ( ncdf_isvar(TRIM(tdir)//'ol_data_used') ) THEN
        CALL ncdf_getvar(TRIM(tdir)//'ol_data_used', readbyte1) 
...

by

      IF ( ncdf_isvar(TRIM(tdir)//'ol_data_available') ) THEN 
        CALL ncdf_getvar(TRIM(tdir)//'ol_data_available', readbyte1)
      ELSE IF (data%leo_id(1:3) == 'SE6') THEN
...

to get a match with Stig's 'before' code. For the record, ol_data_used is present in sx6_test.n4 (the SX-6 dataset attached above) but not in eum_test.n4 (the standard ROPP Metop lev1a test dataset), and ol_data_available is present in eum_test.n4 but not in sx6_test.n4. This is towards the end of the data reading section (Sec 9.8), and as far as I can see it only affects the setting of the PCD_open_loop bit of data%PCD.

[Ian Culverwell]

test_spdp.f90

[Ian Culverwell]

The r6895 mods (the ones to read lev1a data) have been included in r7095. This works OK on the pmode=NRT dataset:

---------------------------------------------------------------------
                  EUMETSAT to ROPP netCDF Converter
---------------------------------------------------------------------

INFO (from eum2ropp):  Reading file ../data/sx6_test.n4
INFO (from ropp_io_read_eumdata):  Format Version 14.0
 pmode = NRT                           
 frequencies =    1604812500.0000000        1248187500.0000000     
INFO (from eum2ropp):  Profile    1 : OC_20200111011620_SE6A_R003_EUME
INFO (from eum2ropp):  Writing sx6_test_l.nc

The printed out frequencies match those in the input file:

idculv@vld189:> ncks -H -Q -v /data/level_1a/open_loop/frequencies ../data/sx6_test.n4
netcdf sx6_test {
  group: data {
    group: level_1a {
      group: open_loop {
        dimensions:
          codes = 2 ;

        variables:
          string codes(codes) ;

          double frequencies(codes) ;

        data:
          codes = "1c", "2c" ;

          frequencies = 1604812500, 1248187500 ;

      } // group /data/level_1a/open_loop
    } // group /data/level_1a
  } // group /data
} // group /

and the output file

idculv@vld189:> ncks -H -Q -vfreq_L1,freq_L2 sx6_test_l.nc
netcdf sx6_test_l {
  dimensions:
    dim_unlim = UNLIMITED ; // (1 currently)

  variables:
    double freq_L1(dim_unlim) ;

    double freq_L2(dim_unlim) ;

  data:
    freq_L1 = 1604812500 ;

    freq_L2 = 1248187500 ;

} // group /

[Ian Culverwell]

S6A_RO_1B_BND_20220201T222900_20220201T223046_20220213T151508_G26NNEUMOPE_NT_R02.nc

[Ian Culverwell]

S6A_RO_1B_BND_20220201T235335_20220201T235516_20220213T152117_R20NNEUMOPE_NT_R02.nc

[Ian Culverwell]

Stig has sent two newer example SX6 files, one for GPS and one for GLONASS. Both attached. These include /data/level_1a/combined/codes:

idculv@vld189:> ncks -H -Q -v/data/level_1a/combined/codes $file
netcdf S6A_RO_1B_BND_____20220201T235335_20220201T235516_20220213T152117_R20__NN__________EUM__OPE_NT_R02 {
  group: data {
    group: level_1a {
      group: combined {
        dimensions:
          codes = 2 ;

        variables:
          string codes(codes) ;

        data:
          codes = "1c", "2c" ;

      } // group /data/level_1a/combined
    } // group /data/level_1a
  } // group /data
} // group /

It would be good to read these codes. Unfortunately the current ROPP netCDF reader cannot read arrays of strings. Look into this.

[Ian Culverwell]

Both these can be read into ROPP, and the frequencies are correctly ingested:

          frequency_1c = 1575420000 ;
          frequency_2l = 1227600000 ;

from the GPS file (G20), and

          frequency_1c = 1603125000 ;
          frequency_2c = 1246875000 ;

for the GLONASS file (R20).

[Ian Culverwell]

Struggling to work out how to read the string array codes. We can't even compile the code

       CALL ncdf_getvar('/data/level_1a/combined/codes', codes)

because

gfortran -I../build -O2 -I/data/users/idculv/ROPP/gfortran/include -c -o ropp_io_read_ncdf_get.o  `test -f '../ropp/ropp_io_read_ncdf_get.f90' || echo './'`../ropp/ropp_io_read_ncdf_get.f90
../ropp/ropp_io_read_ncdf_get.f90:3209.64:

        CALL ncdf_getvar('/data/level_1a/combined/codes', codes)
                                                                1
Error: There is no specific subroutine for the generic 'ncdf_getvar' at (1)
make: *** [ropp_io_read_ncdf_get.o] Error 1

I think we would need to write a new routine ncdf_putgetvar_achar.m4, to read arrays of strings. I might just have managed this, if we hadn't been forced to muck around with ncdf_putgetvar_schar.m4 (the corresponding code for scalar strings) at ROPP-10.0, so that it could cope with strings rather than characters. This meant introducing C code from Leonid at EUMETSAT, which is hard enough to understand anyway, let alone generalise to arrays of strings.

Speak to Chris M about this.

[Ian Culverwell]

There is one possible workaround. The /data/occultation/gnss_system variable says which GNSS satellite transmitted the message, e.g.

 gnss_system = Glonass

for the Glonass (R20) file, and

 gnss_system = GPS

for the Metop (G26) file. Is this enough to tell us the codes that will be in the file?

[Ian Culverwell]

Speak to Stig about this

[Ian Culverwell]

Work on the other bits of r6738.

gfz2ropp

Add GPS freqs to ROPP output:

  ROdata%Signal1%Freq = 1575420000d0  ! L1 GPS frequency
  ROdata%Signal2%Freq = 1227600000d0  ! L2 GPS frequency

These are output OK:

 freq_L1 = 1.5754e+09 ;
 freq_L2 = 1.2276e+09 ;

Committed change at r7096.

[Ian Culverwell]

ropp_fm_bg2ro_1d

The signal frequencies are copied from the levels file in ropp_fm_bg2ro_1d, if requested:

  ro_data%signal1 = l_data%signal1
  ro_data%signal2 = l_data%signal2

Builds OK, and all FM core tests pass, so commit at r7097.

[Ian Culverwell]

Asked CM about reading vectors of strings from netCDF files. He'll ask Leonid to look at it. He thinks it might take a few days. Meanwhile, using gnss_system to infer the codes should work OK for SX6 data, but it might go wrong for Spire, apparently.

[Ian Culverwell]

For the record: initial core tests of these changes cause segmentation faults when running ropp_pp_invert_tool and ropp_pp_grasrs2ropp.

The ropp_pp_invert_tool fails with:

Testing ROPP PP invert tool...

 
Running t_pp_invert_1 (PP invert; default options) ...

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0x7FA3D542D6D7
#1  0x7FA3D542DD1E
#2  0x7FA3D49223FF
#3  0x40BB4D in ropp_pp_tdry_
#4  0x406C41 in MAIN__ at ropp_pp_invert_tool.f90:?
./test_pp_invert.sh: line 32: 18936 Segmentation fault      ./$EXEC ${EXTRA_CMD} $IFILE -o $OFILE -d >> $LOGFILE
make: *** [test_pp] Error 1

It turns out that this is because

 out_refrac%alt_refrac =                        NaN                       NaN ...
 out_refrac%refrac =                        NaN                       NaN ...

which (understandably) give problems to ropp_pp_tdry.f90. The NaNs in turn are caused by

 ro_data%signal1%freq =   -99999000.000000000     
 ro_data%signal2%freq =   -99999000.000000000     

which causes problems (I infer) in the now frequency-dependent ropp_pp_ionospheric_correction.f90, which calculates out_ba%impact_opt and out_ba%bangle_opt, which are then used to calculate out_refrac%refrac and out_refrac%alt_refrac. We therefore need to default the signal frequencies (to GPS L1 and L2), when reading 1D and 2D ROPP-format files in ropp_io_read_ncdf_get.f90, by replacing

IF ( ncdf_isvar('freq_L1') ) &
  CALL ncdf_getvar('freq_L1', data%signal1%freq,    &
                   units = data%signal1%units%freq, &
                   range = data%signal1%range%freq, &
                   rec   = irec)

(+sim. for L2) by

IF ( ncdf_isvar('freq_L1') ) THEN
  CALL ncdf_getvar('freq_L1', data%signal1%freq,    &
                   units = data%signal1%units%freq, &
                   range = data%signal1%range%freq, &
                   rec   = irec)
ELSE
  data%signal1%freq = 1.57542E9_wp ! GPS L1
ENDIF

(+sim. for L2).

This results in

 ro_data%signal1%freq =    1575420000.0000000     
 ro_data%signal2%freq =    1227600000.0000000     
 out_refrac%alt_refrac =    927.51184770651162        1061.8983530187979 ...
 out_refrac%refrac =    324.92429264108227        319.50986306528148 ...

and means that ROPP passes the ropp_pp_invert_tool core test.

[Ian Culverwell]

It's a similar story for ropp_pp_grasrs2ropp.f90, which fails with

Testing ROPP PP raw sampling tool...

 
Running t_pp_rs_1 (PP raw sampling; default options) ...
*** Error in `./../tools/ropp_pp_occ_tool': free(): invalid next size (fast): 0x0000000001c26550 ***
======= Backtrace: =========
/usr/lib64/libc.so.6(+0x81329)[0x7fb5d0072329]
./../tools/ropp_pp_occ_tool[0x46d6f6]
./../tools/ropp_pp_occ_tool[0x46ed81]
./../tools/ropp_pp_occ_tool[0x45a06a]
./../tools/ropp_pp_occ_tool[0x411f73]
./../tools/ropp_pp_occ_tool[0x40681f]
./../tools/ropp_pp_occ_tool[0x40395d]
/usr/lib64/libc.so.6(__libc_start_main+0xf5)[0x7fb5d0013555]
./../tools/ropp_pp_occ_tool[0x403986]
======= Memory map: ========
00400000-00ace000 r-xp 00000000 00:32 122398829521                       /net/data/users/idculv/ROPP/ropp_src/branches/dev/Share/ROPP111_SX6_read/ropp_pp/tools/ropp_pp_occ_tool
00ccd000-00cce000 r--p 006cd000 00:32 122398829521                       /net/data/users/idculv/ROPP/ropp_src/branches/dev/Share/ROPP111_SX6_read/ropp_pp/tools/ropp_pp_occ_tool
00cce000-00ce1000 rw-p 006ce000 00:32 122398829521                       /net/data/users/idculv/ROPP/ropp_src/branches/dev/Share/ROPP111_SX6_read/ropp_pp/tools/ropp_pp_occ_tool
00ce1000-00f6c000 rw-p 00000000 00:00 0 
01a50000-0415d000 rw-p 00000000 00:00 0                                  [heap]
7fb5c8000000-7fb5c8021000 rw-p 00000000 00:00 0 
7fb5c8021000-7fb5cc000000 ---p 00000000 00:00 0 
7fb5cfff1000-7fb5d01b5000 r-xp 00000000 fd:00 8810066                    /usr/lib64/libc-2.17.so
7fb5d01b5000-7fb5d03b4000 ---p 001c4000 fd:00 8810066                    /usr/lib64/libc-2.17.so
7fb5d03b4000-7fb5d03b8000 r--p 001c3000 fd:00 8810066                    /usr/lib64/libc-2.17.so
7fb5d03b8000-7fb5d03ba000 rw-p 001c7000 fd:00 8810066                    /usr/lib64/libc-2.17.so
7fb5d03ba000-7fb5d03bf000 rw-p 00000000 00:00 0 
7fb5d03bf000-7fb5d03fe000 r-xp 00000000 00:31 13835060284415658851       /net/project/ukmo/rhel7/fortran/opt/gfortran/packages/gcc/8.1.0/lib64/libquadmath.so.0.0.0
7fb5d03fe000-7fb5d05fd000 ---p 0003f000 00:31 13835060284415658851       /net/project/ukmo/rhel7/fortran/opt/gfortran/packages/gcc/8.1.0/lib64/libquadmath.so.0.0.0
7fb5d05fd000-7fb5d05fe000 r--p 0003e000 00:31 13835060284415658851       /net/project/ukmo/rhel7/fortran/opt/gfortran/packages/gcc/8.1.0/lib64/libquadmath.so.0.0.0
7fb5d05fe000-7fb5d05ff000 rw-p 0003f000 00:31 13835060284415658851       /net/project/ukmo/rhel7/fortran/opt/gfortran/packages/gcc/8.1.0/lib64/libquadmath.so.0.0.0
7fb5d05ff000-7fb5d0616000 r-xp 00000000 00:31 13835060284408686368       /net/project/ukmo/rhel7/fortran/opt/gfortran/packages/gcc/8.1.0/lib64/libgcc_s.so.1
7fb5d0616000-7fb5d0815000 ---p 00017000 00:31 13835060284408686368       /net/project/ukmo/rhel7/fortran/opt/gfortran/packages/gcc/8.1.0/lib64/libgcc_s.so.1
7fb5d0815000-7fb5d0816000 r--p 00016000 00:31 13835060284408686368       /net/project/ukmo/rhel7/fortran/opt/gfortran/packages/gcc/8.1.0/lib64/libgcc_s.so.1
7fb5d0816000-7fb5d0817000 rw-p 00017000 00:31 13835060284408686368       /net/project/ukmo/rhel7/fortran/opt/gfortran/packages/gcc/8.1.0/lib64/libgcc_s.so.1
7fb5d0817000-7fb5d0918000 r-xp 00000000 fd:00 8862654                    /usr/lib64/libm-2.17.so
7fb5d0918000-7fb5d0b17000 ---p 00101000 fd:00 8862654                    /usr/lib64/libm-2.17.so
7fb5d0b17000-7fb5d0b18000 r--p 00100000 fd:00 8862654                    /usr/lib64/libm-2.17.so
7fb5d0b18000-7fb5d0b19000 rw-p 00101000 fd:00 8862654                    /usr/lib64/libm-2.17.so
7fb5d0b19000-7fb5d0c38000 r-xp 00000000 fd:00 8389555                    /usr/lib64/libgfortran.so.3.0.0
7fb5d0c38000-7fb5d0e38000 ---p 0011f000 fd:00 8389555                    /usr/lib64/libgfortran.so.3.0.0
7fb5d0e38000-7fb5d0e39000 r--p 0011f000 fd:00 8389555                    /usr/lib64/libgfortran.so.3.0.0
7fb5d0e39000-7fb5d0e3b000 rw-p 00120000 fd:00 8389555                    /usr/lib64/libgfortran.so.3.0.0
7fb5d0e3b000-7fb5d0e3d000 r-xp 00000000 fd:00 8845516                    /usr/lib64/libdl-2.17.so
7fb5d0e3d000-7fb5d103d000 ---p 00002000 fd:00 8845516                    /usr/lib64/libdl-2.17.so
7fb5d103d000-7fb5d103e000 r--p 00002000 fd:00 8845516                    /usr/lib64/libdl-2.17.so
7fb5d103e000-7fb5d103f000 rw-p 00003000 fd:00 8845516                    /usr/lib64/libdl-2.17.so
7fb5d103f000-7fb5d1061000 r-xp 00000000 fd:00 8405752                    /usr/lib64/ld-2.17.so
7fb5d1237000-7fb5d123c000 rw-p 00000000 00:00 0 
7fb5d125d000-7fb5d1260000 rw-p 00000000 00:00 0 
7fb5d1260000-7fb5d1261000 r--p 00021000 fd:00 8405752                    /usr/lib64/ld-2.17.so
7fb5d1261000-7fb5d1262000 rw-p 00022000 fd:00 8405752                    /usr/lib64/ld-2.17.so
7fb5d1262000-7fb5d1263000 rw-p 00000000 00:00 0 
7ffc1b32d000-7ffc1b352000 rw-p 00000000 00:00 0                          [stack]
7ffc1b3ac000-7ffc1b3ae000 r-xp 00000000 00:00 0                          [vdso]
ffffffffff600000-ffffffffff601000 r-xp 00000000 00:00 0                  [vsyscall]

Program received signal SIGABRT: Process abort signal.

Backtrace for this error:
#0  0x7FB5D0B326D7
#1  0x7FB5D0B32D1E
#2  0x7FB5D00273FF
#3  0x7FB5D0027387
#4  0x7FB5D0028A77
#5  0x7FB5D0069F66
#6  0x7FB5D0072328
#7  0x46D6F5 in ropp_pp_radiooptic_analysis_
#8  0x46ED80 in ropp_pp_radioholographic_filter_
#9  0x45A069 in ropp_pp_correct_l2_
#10  0x411F72 in ropp_pp_preprocess_
#11  0x40681E in MAIN__ at ropp_pp_occ_tool.f90:?
./test_pp_rs.sh: line 33: 116675 Aborted                 ./$EXEC2 ${EXTRA_CMD} $TFILE -o $OFILE -c $CFILE -m $MOPT -d >> $LOGFILE
make: *** [test_pp] Error 1

I didn't do much investigation into this, but adding

          ! 6.4 Signal characteristics

          ro_data%signal1%freq = 1.57542E9_wp ! GPS L1
          ro_data%signal2%freq = 1.2276E9_wp  ! GPS L2

to ropp_pp_grasrs2ropp.f90 allows ROPP to pass the ropp_pp_grasrs2ropp core test.

[Ian Culverwell]

Both of the above two changes were included in r7101, and allow the ropp_pp module to pass all its core tests:

************************** SUMMARY OF ROPP_PP TEST RESULTS ***************************
--------------------------------------------------------------------------------------
|                   Test name    |              Description       |    Run? |  PASS? |
--------------------------------------------------------------------------------------
|                  t_pp_invert_1 |     PP invert; default options |     Run |  PASS  |
|                     t_pp_occ_1 |        PP occ; default options |     Run |  PASS  |
|                t_pp_occ_gnos_1 |              PP occ; GNOS data |     Run |  PASS  |
|                      t_pp_rs_1 | PP raw sampling; default optio |     Run |  PASS  |
|                    t_pp_abel_1 |          PP Abel/Inv; def opts |     Run |  PASS  |
|                t_pp_spectra_1a |    PP spectra; def opt (L1 dt) |     Run |  PASS  |
|                t_pp_spectra_1b |    PP spectra; def opt (L2 dt) |     Run |  PASS  |
|                t_pp_spectra_1c |    PP spectra; def opt (L1 ep) |     Run |  PASS  |
|                t_pp_spectra_1d |    PP spectra; def opt (L2 ep) |     Run |  PASS  |
|                    t_pp_wopt_1 |      PP 1D WOPT; quick options |     Run |  PASS  |
|                 t_pp_wopt_2d_1 |      PP 2D WOPT; quick options |     Run |  PASS  |
--------------------------------------------------------------------------------------

[Ian Culverwell]

DMI's Sentinel-6 config file (from their r6831) added to branch at r7102.

[Ian Culverwell]

Pages 22-27 of ​http://acc.igs.org/misc/rinex304.pdf list the possible carrier phases (frequencies and 'modes') for the GNSS transmitters we are interested in.

[Ian Culverwell]

A first go at including residual ("kappa") bending angle correction in ROPP has been made at r7103. There are three issues/problems:

• cin(1) is no longer redefined (by means of cin(1) = cin(1)*2.0_wp*f_L1**4 / (f_L1**2 - f_L2**2)**2) if bangle_kappa is not used. Stig says this is deliberate, because the original rescaling was apparently a bug, and results in less weight being given to the data than the config%so_method = so option. See r5646. This will change default behaviour.
• icdr_mode = .FALSE. by default. This means that the kappa correction is applied. This would change default behaviour. Stig is OK with me changing the default, as long as it's called something other than icdr_mode, which I support. (I was thinking of adding a configuration file option called something like kappa_corr, which would be .FALSE. by default.)
• Need to document this in the ROPP PP User Guide.

[Ian Culverwell]

Following discussions with Stig, we have agreed that it is OK to allow the rescaling of cin(1) if config%so_method = lcso. My take on it is this.

The ropp_pp_ionospheric_correction.f90 lines

  ! 7.3 Estimation of ionospheric noise covariance
 
  ba_neut(i_ltr:nh) = ba_low(1,i_ltr:nh) -   &
                              ba_low(2,i_ltr:nh) * (f_L2/f_L1)**2      
 
  IF (iil > 0) THEN
     cin(1) = (SUM(ba_neut(iil:iiu)**2)/(iiu-iil+1))/2.0_wp
     cin(2) = ((f_L1/f_L2)**4)*(SUM(ba_neut(iil:iiu)**2)/(iiu-iil+1))/2.0_wp
  ELSE
     cin(:) = 0.0_wp
  ENDIF

mean that ba_neut equals (1 – f2²/f1²) Δα_N, where Δα_N is the smoothed neutral bending angle, or at least the bit of it that remains after the model bangle (and the kappa correction) have been subtracted off. We then get cin(1) = (1/2) (1 – f2²/f1²)² (mean square of Δα_N). I guess this is taken as an estimate of the ionospheric noise covariance because the index range (iil:iiu) is in the ionosphere. The formula cin(2) = cin(1) * (f1/f2)⁴ suggests to me that cin(1) is the noise on L1, because we might suspect that the standard deviation of the noise at the two frequencies is inversely proportional to the frequency squared. The lcso rescaling then reduces cin(1) to simply 1*(variance of Δα_N), which sounds plausible to me, because you add it to the (neutral) model bangle_LM multiplied by the neutral signal covariance. Numerically, the cin(1) rescaling factor 2.0_wp*f_L1**4 / (f_L1**2 - f_L2**2)**2 is about 13 (>1), so

    ba_LC(:) = ba_LM(:) + (ba_LC(:)-ba_LM(:))*sem*ba_LM(:)**2 /   &
                  (sem*ba_LM(:)**2 + cin(1))

means that rescaling (i.e. increasing) cin(1) makes ba_LC closer to ba_LM, i.e. it gives less weight to the data.

[Ian Culverwell]

I have also removed the icdr stuff at r7110. (Stig is happy with this.) The kappa correction is now applied if the (new) configuration option kappa_corr = .true.. (It is .false. by default.) These changes were introduced at r7108.

[Ian Culverwell]

ropp_pp_test_occ.png

[Ian Culverwell]

ropp_pp_test_gnos.png

[Ian Culverwell]

ropp_pp_test_rs.png

[Ian Culverwell]

Tried to include the calls to ropp_pp_kappa_residual in ropp_pp_occ_tool.f90. Understandably, this causes 3 of the core tests to fail:

************************** SUMMARY OF ROPP_PP TEST RESULTS ***************************
--------------------------------------------------------------------------------------
|                   Test name    |              Description       |    Run? |  PASS? |
--------------------------------------------------------------------------------------
|                  t_pp_invert_1 |     PP invert; default options |     Run |  PASS  |
|                     t_pp_occ_1 |        PP occ; default options |     Run | *FAIL* |
|                t_pp_occ_gnos_1 |              PP occ; GNOS data |     Run | *FAIL* |
|                      t_pp_rs_1 | PP raw sampling; default optio |     Run | *FAIL* |
|                    t_pp_abel_1 |          PP Abel/Inv; def opts |     Run |  PASS  |
|                t_pp_spectra_1a |    PP spectra; def opt (L1 dt) |     Run |  PASS  |
|                t_pp_spectra_1b |    PP spectra; def opt (L2 dt) |     Run |  PASS  |
|                t_pp_spectra_1c |    PP spectra; def opt (L1 ep) |     Run |  PASS  |
|                t_pp_spectra_1d |    PP spectra; def opt (L2 ep) |     Run |  PASS  |
|                    t_pp_wopt_1 |      PP 1D WOPT; quick options |     Run |  PASS  |
|                 t_pp_wopt_2d_1 |      PP 2D WOPT; quick options |     Run |  PASS  |
--------------------------------------------------------------------------------------

For the record, here are the bangle_kappa corrections that are applied at various stages in the three failing cases:

• t_pp_occ_1

• t_pp_occ_gnos_1

• t_pp_rs_1

These look very 'scratchy' to me, and quite big in places. Mmm.

(Note that to print out the last of these (the raw sampling one), I had to comment out

    DEALLOCATE(ro_data%vlist%VlistD1d%data)

from Sec 2 of ropp_pp_preprocess_grasrs.f90. Without doing this, I couldn't add the new bangle_kappa variables to ro_data in ropp_pp_occ_tool.f90, and then plot them. (I got a segmentation fault.) I also had to comment out

        CALL ropp_io_free(ro_data%vlist)   !! interim, avoid writing lcf data

from Sec 22 of ropp_pp_occ_tool.f90. Someone should probably return to this.)

[Ian Culverwell]

I can get around the above outputting problem, without having to comment out the

  DEALLOCATE(ro_data%vlist%VlistD1d%data)

statement in ropp_pp_preprocess_grasrs.f90 (and its counterpart in Sec 3 of ropp_pp_preprocess_cosmic.f90), by adding the line

  IF (ASSOCIATED(ro_data%vlist%VlistD1d)) CALL ropp_io_free(ro_data%vlist)

in ropp_pp_occ_tool.f90 before I add the first new diagnostic.

[Ian Culverwell]

ropp_pp_test_occ_Sec11.png

[Ian Culverwell]

ropp_pp_test_occ_Sec14.png

[Ian Culverwell]

ropp_pp_test_gnos_Sec11.png

[Ian Culverwell]

ropp_pp_test_gnos_Sec14.png

[Ian Culverwell]

ropp_pp_test_rs_Sec11.png

[Ian Culverwell]

ropp_pp_test_rs_Sec14.png

[Ian Culverwell]

The oddest and biggest bangle_kappas occur in Sec 11 of ropp_pp_occ_tool.f90, where it is added to the LC of the GO-processed L1 and L2 signals. This is quite noisy, as is bangle_kappa, which I suppose is reasonable, given their common formulation in terms of differences between noisy signals. And adding the latter doesn't really help the former.

For example, for the standard occ core test:

• t_pp_occ_1

For the GNOS test:

• t_pp_occ_gnos_1

For the RS test:

• t_pp_occ_rs_1

In summary, the LC of the GO bangles is noisy, and so is the kappa correction, and in the same places - but not enough to reverse the noise in the input. (Not surprising - it's not designed to.)

[Ian Culverwell]

Sec 14 of ropp_pp_occ_tool.f90 does a similar thing with the WO-processed LC of bangle_L1 and bangle_L2. These are much smoother, and the kappa correction is much smaller, so adding it has very little effect on the output of a single test. (It should improve climatological/reanalysis averages, however.)

Plots are attached to this ticket but not displayed.

[Ian Culverwell]

ropp_pp_test_invert_Sec10.png

[Ian Culverwell]

Tried to do the same thing with ropp_pp_invert_tool.f90. I think the only code that needs changing is in Sec 10, which calculates the 'LC' of bangle_L1 and bangle_L2. This requires config%method to be set to 'NONE' (by default it is 'GMSIS'). This cannot have been tried very often, because I found I needed to add the lines

          Imin = SUM(MINLOC(ABS(out_ba%impact)))
          Imax = SUM(MAXLOC(ABS(out_ba%impact)))

to that section, as otherwise those variables are not set, which means that the later (Sec 14) lines

        IF (config%output_lev1b) THEN
          out_ba%npoints = imax - imin + 1
          CALL ropp_io_roprof2roprof(out_ba, ro_data%lev1b)
        ELSE
          CALL ropp_io_free(ro_data%lev1b)
        END IF

        IF (config%output_lev2a) THEN
          out_refrac%npoints = imax - imin + 1
          CALL ropp_io_roprof2roprof(out_refrac, ro_data%lev2a)
        ELSE
          CALL ropp_io_free(ro_data%lev2a)
        END IF

prevent lev1b and lev2a data from being added to the output file. (With gfortran I found that without initialisation these variables took the values

imin, imax =            0  -843944955

and I was lucky to avoid a segmentation fault.)

Anyway, assuming this is the right thing to do (check with Stig), the corresponding bangle_kappa are:

• t_pp_invert_1

It's pretty small (~10^-7 rad), and only has a 'spike' where bangle_L1 has a spike.

I think this looks OK.

[Ian Culverwell]

Changes to ropp_pp_occ_tool.f90 have been committed at r7113.

Changes to ropp_pp_invert_tool.f90 have been committed at r7114.

[Ian Culverwell]

Cosmetic changes to ropp_pp_ionospheric_correction.f90 committed at r7115.

[Ian Culverwell]

ropp_pp_test_invert_Sec11.png

[Ian Culverwell]

Stig thinks (and I now agree) that the kappa correction should also be applied to the statistically optimised bending angles bangle in Sec. 11.1 of ropp_pp_invert_tool.f90. For the standard test file the correction is about the same size as it was in Sec 10, which makes sense:

• t_pp_invert_1

Again, I think this is OK.

[Ian Culverwell]

This change to ropp_pp_invert_tool.f90 was made at r7120.

[Ian Culverwell]

Documentation of kappa correction added to the ROPP-11.1 PP UG at r7137.

[Ian Culverwell]

r7092:r7140 of https://trac.romsaf.org/ropp/browser/ropp_src/branches/dev/Share/ROPP111_SX6_read have been merged into https://trac.romsaf.org/ropp/browser/ropp_src/branches/dev/Share/ROPP111_prototype at r7141. Any further development should take place on that branch, so I'm closing this ticket.