id summary reporter owner description type status priority milestone component version resolution keywords cc 661 Old bug in ropp_fm_state2state_ecmwf_ad.f90 Stig Syndergaard Ian Culverwell "The adjoint (AD) of the interpolation coefficients in ropp_fm_state2state_ecmwf_ad.f90 includes the following line: {{{ p_hlv_ad(1:n_hlv-1) = p_hlv_ad(1:n_hlv-1) - alpha_ad/del_p*ln_prflv }}} The tangent linear (TL) of the interpolation coefficients in ropp_fm_state2state_ecmwf_tl.f90 includes this line: {{{ alpha_tl = del_p_tl*p_hlv(2:n_hlv)/(del_p*del_p)*ln_prflv & - ln_prflv_tl*(p_hlv(2:n_hlv)/del_p) & - p_hlv_tl(2:n_hlv)/del_p*ln_prflv }}} These two are inconsistent. The tangent linear is correct according to the expression in eq. (3.9) in Simmons and Burridge (1981): https://journals.ametsoc.org/doi/pdf/10.1175/1520-0493%281981%29109%3C0758%3AAEAAMC%3E2.0.CO%3B2 The correct line in ropp_fm_state2state_ecmwf_ad.f90 should therefore be: {{{ p_hlv_ad(2:n_hlv) = p_hlv_ad(2:n_hlv) - alpha_ad/del_p*ln_prflv }}} which can be verified by writing matrices on a piece of paper or simply following the rules of how to write AD code based on TL code. A fix has been implemented in the dmi_trunk_9.0 at r6115, although we do not at DMI use the AD code in the 1Dvar (since we are using the LEVMARQ method and the AD code is only used with the MINROPP method). It is unclear what other users could be affected. The bug has been in ROPP from day 1 as far as I can see (from looking at the svn history). If the code originally comes from ECMWF, it would probably be a good idea for some to check if such code there has the same issue. Before the fix (actually unfixing it), the following results were obtained using the test code made by Joe here: https://trac.romsaf.org/ropp/browser/ropp_src/branches/dev/Share/jkn_dmi_trunk_9.0.2/ropp_fm/tests/t_fascod_ecmwf_ad.f90 {{{ AD: x%non_ideal, x%new_ref_op, x%new_bangle_op = F F F Checking REFRACTIVITY AD FASCOD_TRO norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 364528929.066 364528929.045 1.0000000 0.57491E-10 Checking BENDING ANGLE AD FASCOD_TRO norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = -9.025 -9.025 1.0000000 -0.12613E-09 Checking REFRACTIVITY AD FASCOD_MLS norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 325588949.035 325588949.019 1.0000000 0.49311E-10 Checking BENDING ANGLE AD FASCOD_MLS norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 75.279 75.279 1.0000000 0.16198E-09 Checking REFRACTIVITY AD FASCOD_MLW norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 431837535.061 431837535.253 1.0000000 -0.44447E-09 Checking BENDING ANGLE AD FASCOD_MLW norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = -79.704 -79.704 1.0000000 0.59795E-11 Checking REFRACTIVITY AD FASCOD_SAS norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 372208695.410 372208695.450 1.0000000 -0.10889E-09 Checking BENDING ANGLE AD FASCOD_SAS norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 25.611 25.611 1.0000000 0.57009E-10 Checking REFRACTIVITY AD FASCOD_SAW norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 409038225.382 409038225.393 1.0000000 -0.26912E-10 Checking BENDING ANGLE AD FASCOD_SAW norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 8.828 8.828 1.0000000 0.17541E-09 **************************** ********** PASS ********** **************************** }}} With the fix implemented, the test results are: {{{ AD: x%non_ideal, x%new_ref_op, x%new_bangle_op = F F F Checking REFRACTIVITY AD FASCOD_TRO norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 364528929.066 364528929.066 1.0000000 0.16351E-15 Checking BENDING ANGLE AD FASCOD_TRO norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = -9.025 -9.025 1.0000000 -0.84639E-14 Checking REFRACTIVITY AD FASCOD_MLS norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 325588949.035 325588949.035 1.0000000 -0.36613E-15 Checking BENDING ANGLE AD FASCOD_MLS norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 75.279 75.279 1.0000000 0.15102E-14 Checking REFRACTIVITY AD FASCOD_MLW norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 431837535.061 431837535.061 1.0000000 0.00000E+00 Checking BENDING ANGLE AD FASCOD_MLW norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = -79.704 -79.704 1.0000000 -0.30310E-14 Checking REFRACTIVITY AD FASCOD_SAS norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 372208695.410 372208695.410 1.0000000 0.16014E-15 Checking BENDING ANGLE AD FASCOD_SAS norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 25.611 25.611 1.0000000 -0.69360E-15 Checking REFRACTIVITY AD FASCOD_SAW norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 409038225.382 409038225.382 1.0000000 0.00000E+00 Checking BENDING ANGLE AD FASCOD_SAW norm1, norm2, norm1/norm2, (norm1-norm2)/norm2 = 8.828 8.828 1.0000000 0.54330E-14 **************************** ********** PASS ********** **************************** }}} Thus an improvement from about 10^-9^ to 10^-14^. The change does not seem to break anything: doing 'make test' in dmi_trunk_9.0 gives the same results before and after the fix: {{{ ************************** SUMMARY OF ROPP_FM TEST RESULTS *************************** -------------------------------------------------------------------------------------- | Test name | Description | Run? | PASS? | -------------------------------------------------------------------------------------- | t_fascod_1 | FM FASCOD 1D | Run | *FAIL* | | t_fascod_tl_1 | FM_TL FASCOD 1D | Run | PASS | | t_fascod_ad_1 | FM_AD FASCOD 1D | Run | PASS | | t_fascod_2 | FM FASCOD 1D -comp | Run | *FAIL* | | t_fascod_tl_2 | FM_TL FASCOD 1D -comp | Run | PASS | | t_fascod_ad_2 | FM_AD FASCOD 1D -comp | Run | PASS | | t_fascod_3 | FM FASCOD 1D -new_op | Run | *FAIL* | | t_fascod_tl_3 | FM_TL FASCOD 1D -new_op | Run | PASS | | t_fascod_ad_3 | FM_AD FASCOD 1D -new_op | Run | PASS | | t_fascod_4 | FM FASCOD 1D -comp -new_op | Run | *FAIL* | | t_fascod_tl_4 | FM_TL FASCOD 1D -comp -new_op | Run | PASS | | t_fascod_ad_4 | FM_AD FASCOD 1D -comp -new_op | Run | PASS | | t_twodop_1 | FM TWOD | Run | PASS | | t_twodtl_1 | FM_TL TWOD | Run | PASS | | t_twodad_1 | FM_AD TWOD | Run | PASS | | t_twodop_2 | FM TWOD -comp | Run | PASS | | t_twodtl_2 | FM_TL TWOD -comp | Run | PASS | | t_twodad_2 | FM_AD TWOD -comp | Run | PASS | | t_iono_1 | FM L1 and L2 | Run | *FAIL* | | t_iono_tl_1 | FM_TL L1 and L2 | Run | PASS | | t_iono_ad_1 | FM_AD L1 and L2 | Run | PASS | | t_fm_1D_1 | FM 1D; default options | Run | *FAIL* | | t_fm_1D_2 | FM 1D; compress factors | Run | *FAIL* | | t_fm_2D_1 | FM 2D; default options | Run | PASS | | t_fm_2D_2 | FM 2D; compress factors | Run | PASS | | t_fm_iono_1 | FM iono; L_neutral | Run | *FAIL* | | t_fm_iono_2 | FM iono; L1 and L2 | Run | *FAIL* | -------------------------------------------------------------------------------------- }}} The failures here might be because we in the dmi_trunk generally have different codes than in the official ROPP. This will often be the case, and 'make test' in the dmi_trunk cannot generally be expected to pass all checks when comparisons are made to test files from the official ROPP. The failures in the summary above have therefore not been investigated. " defect assigned normal ROPP possible ropp_fm 9.0 ssy@…