I am making my own finite difference computational magnetohydrodynamic code in Fortran 90. Looking at other codes they appear to calculate for example their $x$-derivatives, bb of their variables, e.g. aa like the following:
CALL system_clock(tstart, count_rate)
DO ix = 1, nx - 1
ixp = ix + 1
DO iy = 1, ny
DO iz = 1, nz
bb(ix,iy,iz) = (aa(ixp,iy,iz) - aa(ix,iy,iz)) / dx
END DO
END DO
END DO
CALL system_clock(tstop, count_rate)
PRINT*, "Computation time =", REAL(tstop - tstart) / REAL(count_rate)
In other words they calculate the derivative element by element. Another way of calculating the derivative is array by array, i.e.:
CALL system_clock(tstart, count_rate)
bb(1:nx-1,:,:) = (aa(2:nx,:,:) - aa(1:nx-1,:,:)) / dx
CALL system_clock(tstop, count_rate)
PRINT*, "Computation time =", REAL(tstop - tstart) / REAL(count_rate)
The second method seems to be about 37 times quicker. So why do I not see more codes calculating their $x$-derivatives array by array instead of element by element? What is the downside of array by array, is there something I am missing?