1dvar electron density retrieval code

[Ian Culverwell]

At their meeting on 11/1/2021 the ROPP Coordination Group agreed to defer IEEC’s AVHIROv2/GRID2EPD-based electron density retrieval method to CDOP4, and to implement Sean's variational retrieval method in CDOP3 instead. This decision was endorsed by the Steering Group at their meeting in March 2021 (and activities were replanned accordingly). This will be the major scientific component of ROPP-11.0.

[Ian Culverwell]

Summary

Sean's original code is in the attached sean.tar.gz. compile_slant_tec shows how to build it, and slant_all_sean.py plots the results.

Development work is done in the ROPP-11 branch ​https://svn.romsaf.org/ropp//ropp_src/branches/dev/Share/ROPP110_ne_1dvar.

Aim: to reproduce Sean's results in a branch of ROPP, both for data from IEEC and from the Metop-A extension experiment at EUMETSAT, with comparable timings. This was eventually achieved, as shown in the attached Ne_1dvar.pdf.

As a result, this branch was merged into the ROPP110_prototype branch at r6696, with some further revisions at r6698 and r6699.

Details to follow.

[Ian Culverwell]

Details

1. There are changes to utils, io, fm and 1dvar modules.

2. New ROprof substructure lev2e to hold iono data: {ne, ne_sigma, r_iono, ne_max, ne_max_sigma, r_peak_sigma, h_zero_sigma and h_grad_sigma}. Note that these have different dimensions: the first 3 are on npoints (~100); the last 8 are on nlayers (~1).

3. New state0DFM structure to hold {ne_max, r_peak, h_zero and h_grad} on nlayers.

4. New forward model (and TL) for dbangle = L2-L1, used in 1dvar.

5. New code to calculate ne(r) (and ne_tl) and VTEC from the state vector.

6. New tool, ropp_1dvar_dbangle, written to carry out retrievals. Based on ropp_1dvar_bangle, and uses similar input files.

7. Updated configuration files to include f1, f2, r_gns, r_leo etc. (f1 and f2 not actually used yet.)

8. New 1dvar LM code (only), which calls the new dbangle FM. This has been written to match Sean's minimization code, which means there are some differences to the bangle and refrac minimisers:

1. Include preconditioning (because ne_max ~ 1e12, h_grad ~ 1e-1 ==> ratio of 1e26 in cov matrix).

b.

! Without this, unused observations affect d2J_dx2.
! FIXME: Why don't we (need to) do this for bangle and refrac?
    DO i = 1, nstate
      K(:, i) = K(:, i) * y%weights(:)
    ENDDO

c.

! Use the old, 'standard' (lambda/10) value of dJ_dx on the 
first pass through the LM code.
! This emulates Sean's iono_1dvar.f90 code.
! FIXME: Why don't we do this for bangle or refrac?
    IF ( lmarq ) THEN
      dJ_dx = dJ_dx_old
    ELSE
      dJ_dx   = MATMUL(KO, delta_y) + MATMUL(Bm1, delta_x)
    ENDIF

d.

! Save x%state so that the printed output records a better max(relative change in state).
! FIXME: Why don't we do this for bangle or refrac?
    state_last = x%state

e.

! FIXME: Why don't we do this for bangle and refrac?
! (We currently multiply lambda by 10 in those codes.)
      lambda = 100.0_wp * lambda

f.

! Use the old, large, LM value of lambda on first pass through the 'standard' lambda/10 loop.
! FIXME: Why don't we do this for bangle or refrac?
      IF ( .NOT. lmarq ) lambda = 0.1_wp *lambda

9. Someone who gives a toss might like to examine these differences at some point.

10. New bg and obs data and covariance files, and configuration files.

11. New QC and diagnostic routines, to calculate O-B, A, etc.

12. New I/O routines to read and write the new data.

13. New I/O tools to process ascii IEEC data into ROPP format, and another option (-q) added to eum2ropp to convert Metop-A extension data into ROPP format files.

14. New core test (based on a retrieval using IEEC data) written to automatically test these routines.

15. New code to calculate Gamma (needed for VTEC). (RSR?).

[Ian Culverwell]

Future plans

• Read frequencies from input files.

• Allow user to do retrievals on L1 or L2 as well as L2-L1. (Formally just involves changing the prop. const. in bangle \propto dSTEC/da.)

• Maybe a tool to generate dbangle from state vectors.

[Ian Culverwell]

Further revisions to this code have been made at r6706, r6712, r6715, r6716, r6717, r6731, r6732, r6733, r6736, r6739, r6740, r6741, r6742, r6746, r6751, r6763, r6764 and r6767, mainly as a result of putting ROPP11 through the test folder. (Fortunately the ROM SAF will be getting rid of that soon.)

Most of these are small admin changes but r6739 is a real doozy: it changes

 	  ALLOCATE(Rm1(m, m), Bm1(n, n), K(m, n), K_masked(n, m), KtRm1K(n, n), S(n,n))

to

 	  ALLOCATE(Rm1(m, m), Bm1(n, n), K(m, n), K_masked(m, n), KtRm1K(n, n), S(n,n))

(i.e. it transposes K_masked).

You might have thought that this change was essential, given that the code later says

  K_masked = K

and

  CALL matrix_toast(Rm1, TRANSPOSE(K_masked), KtRm1K)

(i.e. it calculates TRANSPOSE(K_masked) # Rm1 # K_masked) but apparently not: the gfortran, ifort17, nagfor61 and ftn compilers silently reshape K_masked after

  K_masked = K

so that it has the same shape as K, i.e. (m, n). Then everything is fine. Only g95, ifort16 and pgf95 have a problem. g95 crashes out with

Fortran runtime error: Matrix size mismatch in MATMUL()

(which is the only reason I spotted this mistake). ifort16 and pgf95 don't reshape, but still manage to do the matrix toasting, and produce the same results as the reshaping compilers do - presumably the later multiplication of K_masked by obs%weights, and perhaps the fact that m > n, mean that this madness has no effect.

The fact that ifort16 and ifort17 behave differently in this regard is weirdness²!

The attached test program test14.f90 examines this issue. Its output is test14.out.

[Ian Culverwell]

test14.f90

[Ian Culverwell]

test14.out

[Ian Culverwell]

Part of the reason for this is that K and k_masked are ALLOCATABLE. If you dimension

  INTEGER,  PARAMETER                   :: m=3, n=2
  REAL(wp)                              :: K1(n, m), K2(m, n)

in the test program, gfortran (the only one I've tried) kicks up these compile-time errors:

test14.f90:36.2:

  K2 = K1  ! Let's see if this works
  1
Error: Different shape for array assignment at (1) on dimension 
1 (3 and 2)
test14.f90:42.24:

          SHAPE( MATMUL(K1, TRANSPOSE(K2)) )! not defined if Ks 
keep their orig
                        1
Error: Different shape on dimension 2 for argument 'matrix_a' 
and dimension 1 for argument 'matrix_b' at (1) for intrinsic 
matmul
test14.f90:44.49:

  print*, 'MATMUL(K1, TRANSPOSE(K2)) = ', MATMUL(K1, TRANSPOSE(K2))
                                                 1
Error: Different shape on dimension 2 for argument 'matrix_a' 
and dimension 1 for argument 'matrix_b' at (1) for intrinsic 
matmul

as expected.

[Ian Culverwell]

The workaround of r6676 has been reversed at r6902. (Also removed an unnecessary print statement.)

[Ian Culverwell]

This has been incorporated in ROPP-11.0. (In fact, it's the major scientific component.) Passed the DRR, so closing ticket as "done".