Potential bug in sedshi.F90

[jmaerz]

Hi @JorgSchwinger and @TomasTorsvik , while implementing output for burial fluxes, I stumbled over the potential issue described in #27 (not urgent), but in addition in sedshi.F90: could it be that the downward shifting into the burial layer is not correct? To my understanding, the index k should be replaced by ks (loop in l. 135-148 in current master) - note the closure of k-loop before:

            uebers=wsed(i,j)*sedlay(i,j,k,iv)
            sedlay(i,j,ks ,iv)=sedlay(i,j,ks ,iv)-uebers
            burial(i,j,iv)=burial(i,j,iv)+uebers*seddw(k)*porsol(i,j,k)

should read:

            uebers=wsed(i,j)*sedlay(i,j,ks,iv)
            sedlay(i,j,ks ,iv)=sedlay(i,j,ks ,iv)-uebers
            burial(i,j,iv)=burial(i,j,iv)+uebers*seddw(ks)*porsol(i,j,ks)

If so, I would implement the bug-fix it in a branch for beyond-CMIP6 (maybe it would be good to update that branch with the current master before).

[TomasTorsvik]

Hi @jmaerz , I agree that this looks suspicious. Looking at the code history:

https://github.com/NorESMhub/BLOM/blame/master/hamocc/sedshi.F90

most of the sedshi.F90 code has not been changed for a long time, except l. 143 that was included in PR #189.

As it happens, k should take the value ks in this case, since the previous loop terminates at ks-1, but I don't think there is any good reason to rely on this implied behavior.

See discussion in https://stackoverflow.com/questions/62229922/value-of-index-variable-after-loop-exit-in-fortran
and example:

  program main
  implicit none
  integer :: i, n
  n = 10

  do i = 1, n - 1
     print *, "inside: ", i
  enddo

  print *, "after exit: ", i   !! guaranteed to be n?
  end

[TomasTorsvik]

Also, it would be a good idea to update beyond-CMIP6 from master before including the fix.

[JorgSchwinger]

Hi @jmaerz and @TomasTorsvik

It looks like k=ks at this point of the code, right?

In this case, I don't think this classifies as a bug (since it relies on behavior that is part of the Fortran standard). But it is definitely not nice and I agree k should be changed to ks here. There should be no effect on the results (hopefully).

[jmaerz]

Hi @JorgSchwinger an @TomasTorsvik , ok, yes, I see what you mean. Being explicit about it wouldn't harm, though. I guess, in this case, I will add the change simply to my nitrogen cycle branch which will enter the beyone-CMIP6 branch at some point (or the master, if is has been merged).

[TomasTorsvik]

Hi @JorgSchwinger an @TomasTorsvik , ok, yes, I see what you mean. Being explicit about it wouldn't harm, though. I guess, in this case, I will add the change simply to my nitrogen cycle branch which will enter the beyone-CMIP6 branch at some point (or the master, if is has been merged).

Seem like a good plan to me :)