Possible bug in ocn_tracer_advection_std

[cbegeman]

I think that high_order_vert_flux should be allocated to size nVertLevels+1 here so that in the following lines, high_order_vert_flux(k+1) has a reference when maxLevelCell = nVertLevels. @mattdturner Do you agree?

MPAS-Model/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F

Lines 222 to 225 in 29a819a

	do k = 1,maxLevelCell(iCell)
	tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + verticalDivergenceFactor(k) * (high_order_vert_flux(k+1, iCell) &
	- high_order_vert_flux(k, iCell))
	end do ! k loop

[mattdturner]

@cbegeman I agree that this is a bug as the code is currently written. It doesn't look like this particular subroutine is tested in the nightly regression suite (all tests have config_monotonic = .true.), so that would explain why it hasn't shown up as a FAIL before.

I'm not sure, though, if just allocating high_order_vert_flux to size nVertLevels+1 is the correct fix, or if the calculation of tend needs to be modified to handle k == maxLevelCell(iCell) a different way. If we increase the size of high_order_vert_flux by 1, then at the calculation of tend for k == maxLevelCell(iCell) high_order_vert_flux(k+1,iCell) will just have a value of 0.0 (which could very well be the desired behavior).

Perhaps @mark-petersen might have a better understanding of the calculations in this routine, and whether increasing the size of high_order_vert_flux to nVertLevels+1 is the correct resolution to this bug.

[cbegeman]

@mattdturner Thanks for taking a look. My understanding is that high_order_vert_flux(maxLevelCell(iCell)+1,iCell) corresponds to the tracer flux from the seafloor into the bottom cell. I do believe that for most applications this will be equal to 0 but perhaps for some tracers there would be entrainment (Fe, sediment). I'm not sure if that entrainment would be handled here or elsewhere.

[philipwjones]

@mattdturner, @cbegeman Before I modified the mono advection routine, the vertical advection here was essentially identical, so I had planned to make the same performance mods in the std routine that I made in mono. I've been working on a version with GPU mods for both and will be making the changes as part of that, but if we want to make the changes sooner, I can speed that up or someone can copy how the vertical advection is done in the mono routine.

[cbegeman]

@philipwjones No rush on my end. Just wanted to point it out when I came across it. Thanks!