Bug in ropp_fm_interpol.f90

[Ian Culverwell]

Extrapolation below the least input point in ropp_fm_interpol.f90 is incorrect.

  DO k = 1, SIZE(newx)
     j = 1
     DO WHILE (j < SIZE(x) .AND. x(j) < newx(k))
        j = j + 1
     ENDDO
     i = j - 1

     interp(k) = array(j) + &
          ( (newx(k) - x(j)) / (x(i) - x(j)) * (array(i) - array(j)) )
 if (k == 1) print*,'LIN: i, j, array(i), array(j), newx(k), x(i), x(j), array(i), array(j), interp(k) = ', &
                          i, j, array(i), array(j), newx(k), x(i), x(j), array(i), array(j), interp(k)
  ENDDO

generates

 LIN: i, j, array(i), array(j), newx(k), x(i), x(j), array(i), array(j), interp(
 k) =            0           1  4.940656458412465E-324   224.465585702694     
   663.892771230079       0.000000000000000E+000   3684.59747913574     
  4.940656458412465E-324   224.465585702694        40.4443309158695     

Note x(0) = 0.000000000000000E+000 and y(0) = 4.940656458412465E-324. Very dodgy!

It should say, as ropp_fm_interpol_log.f90 does,

  DO k = 1, SIZE(newx)
     j = 2
     DO WHILE (j < SIZE(x) .AND. x(j) < newx(k))
        j = j + 1
     ENDDO
     i = j - 1

     interp(k) = array(j) + &
          ( (newx(k) - x(j)) / (x(i) - x(j)) * (array(i) - array(j)) )
 if (k == 1) print*,'LIN: i, j, array(i), array(j), newx(k), x(i), x(j), array(i), array(j), interp(k) = ', &
                          i, j, array(i), array(j), newx(k), x(i), x(j), array(i), array(j), interp(k)
  ENDDO

This correctly generates

 LIN: i, j, array(i), array(j), newx(k), x(i), x(j), array(i), array(j), interp(
 k) =            1           2   224.465585702694        213.429474246952     
   663.892771230079        3684.59747913574        3702.16006527328     
   224.465585702694        213.429474246952        2122.63898512385     

in this case.

(This has not caused a problem before because only ropp_fm_interpol_log is called in practice.)

[Ian Culverwell]

Note that ropp_fm_interpol_tl.f90 and ropp_fm_interpol_ad.f90 will also need to be updated.

[Ian Culverwell]

Test this by changing

  CALL ropp_fm_interpol_log(z_geop, y%geop, refrac, y%refrac)

to

  CALL ropp_fm_interpol(z_geop, y%geop, refrac, y%refrac)

in ropp_fm/refrac_1d/ropp_fm_refrac_1d.f90, and ditto for _tl and _ad versions. Then call ropp_fm/tests/t_fascod_1.sh (ditto for _tl and _ad). We find that t_fascod_1 fails, as it should, because it's comparing against reference values calculated with log interpolation, but that the _tl and _ad codes test out OK, because they generate the reference data on the fly. (The numbers in the comparison tests are different - correctly - but the threshold criteria are still met.)

Changes committed at r4754.

[Ian Culverwell]

(The test routines t_fascod.f90 (etc) could do with tidying up, but they work and they're only used once, when testing the installation, so leave them be.)

Closing ticket.