This weekend I ran into a rather unexpected issue and I hope you can point out my error. I simply fail to see it and that annoys me.

The program is fairly straightforward: an attempt to solve an advection-diffusion problem with a compact implementation. In the absence of advection, the result should be symmetric.It is if I use method 2 in the code. But with method 1, the result is not and I fail to see why.

Here is the code:

```
! advdiff.f90 --
! Module implementing a straightforeward solver for an advection-diffusion
! equation in two dimensions
!
module advdiff
use iso_c_binding
implicit none
!
! Type holding the various parameters
! (alternatives: put them in separate arguments
! store them in module variables)
!
type, bind(C) :: param_def
real(kind=c_double) :: diffusion ! Diffusion coefficient (m2/s)
real(kind=c_double) :: flow_velocity ! Flow velocity (left to right; m/s)
real(kind=c_double) :: deltax ! Grid cell size (m)
real(kind=c_double) :: deltat ! Time step (s)
integer :: nreport ! Number of steps to set before returning
end type param_def
contains
subroutine update( n, conc_in, conc_out, params ) bind(c, name = 'update')
integer, intent(in) :: n
real(kind=c_double), intent(in) :: conc_in(n,n)
real(kind=c_double), intent(out) :: conc_out(n,n)
type(param_def), intent(in) :: params
integer :: i
real(kind=c_double) :: diffdt, veldt
real(kind=c_double) :: tmp(n-2,n-2)
!
! We want a compact calculation
!
conc_out = conc_in
diffdt = params%diffusion / params%deltax**2 * params%deltat
veldt = params%flow_velocity / params%deltax * params%deltat
!
! Using associate we can formulate the difference equation in a way that avoids
! explicit loops and boundary conditions.
! Note:
! It is not necessarily very efficient ;)
!
! By first copying the initial values and then using the output array, we
! avoid intermediate copying.
!
associate( ccentre => conc_out(2:n-1,2:n-1), &
cleft => conc_out(1:n-2,2:n-1), &
cright => conc_out(3:n ,2:n-1) )
do i = 1,params%nreport
!
! Method 1: produces asymmetric result
!
!!ccentre = ccentre + veldt * ( cleft - ccentre ) + &
!! diffdt * ( cleft + cright - 2.0_c_double * ccentre )
!
! Method 2: produces the expected symmetric result
!
tmp = veldt * ( cleft - ccentre ) + &
diffdt * ( cleft + cright - 2.0_c_double * ccentre )
ccentre = ccentre + tmp
!write(20,'(10f8.4)') tmp
write(20,'(10f8.4)') ccentre
enddo
end associate
end subroutine update
end module advdiff
! test_advdiff.f90 --
! Test program for the advection-diffusion solver
!
program test_advdiff
use advdiff
implicit none
type(param_def) :: params
integer, parameter :: n = 12
integer :: i, j
real(kind=c_double) :: conc_in(n,n), conc_out(n,n)
open( 20, file = 'test_advdiff.out' )
conc_in = 0.0_c_double
conc_in(2:n-1,2:n-1) = 1.0_c_double
params%diffusion = 0.1_c_double
params%flow_velocity = 0.0_c_double
params%deltax = 0.1_c_double
params%deltat = 0.02_c_double
params%nreport = 1
do i = 1,10
call update( n, conc_in, conc_out, params )
write( 20, '(a,i0)' ) 'Time: ', i
do j = 1,n
write( 20, '(100f8.4)' ) conc_out(:,j)
enddo
conc_in = conc_out
enddo
end program test_advdiff
```

The asymmetry occurs right at the first time step.

Also: the results are the same with gfortran and Intel Fortran oneAPI, so I suspect some problem in my code, but I do not see it