Ok,
So here's what being done in the real (obviously illegal code). We are doing advection in x-,y- and z-direction. The advection driver routines all have the structure...
subroutine adv_x(phinew,phiold)
implicit none
real, intent(in) :: phiold(:,:,:)
real, intent(out) :: phinew(:,:,:)
loop over domain
tmp = f(phiold)
loop over domain
phinew = g(tmp)
end subroutine adv_x
...where f() and g() are some complicated explicitly coded out function which is not shown for brevity. The same is true for y- and z-direction (adv_y, adv_z). Note that phi is a rank-3 array. As you can see, the routine is "by design" safe if phiold and phinew correspond to the same memory location. There are two driver routines which call the advection operators, which are called alternatively every timestep with phiold containing the field valid at the last timestep...
subroutine adv_xyz(phinew, phiold)
implicit none
real, intent(in) :: phiold(:,:,:)
real, intent(out) :: phinew(:,:,:)
call adv_x(phinew, phiold)
call adv_y(phinew, phinew)
call adv_z(phinew, phinew)
end subroutine adv_xyz
subroutine adv_zyx(phinew, phiold)
implicit none
real, intent(in) :: phiold(:,:,:)
real, intent(out) :: phinew(:,:,:)
call adv_z(phinew, phiold)
call adv_y(phinew, phinew)
call adv_x(phinew, phinew)
end subroutine adv_xyz
...these could be written "legally", by using...
subroutine adv_xyz(phinew, phiold)
implicit none
real, intent(in) :: phiold(:,:,:)
real, intent(out) :: phinew(:,:,:)
call adv_x(phinew, phiold)
phiold = phinew
call adv_y(phinew, phiold)
phiold = phinew
call adv_z(phinew, phiold)
end subroutine adv_xyz
...but that would degrade performance significantly since this is a memory bound code. To my understanding, the reason that "call adv_x(phinew, (phinew)) works, is that the NAG compiler is generating an extra copy of phinew in memory before calling adv_x, and that is exactly what I would like to avoid. I see the solution of having two implementations of sub_x, -y, -z with one or two arguments, but if there is a significant amount of code inside these subroutines and I will simply have to copy it and change only the interface, which does not help maintainability of the code. So say I declare phinow and phinew as POINTERs inside adv_x, -y, -z. Would that be ok? Would I have to declare them as pointers up through the full calling tree (i.e. also in adv_xyz, -zyx and their callers)?
Cheers,
Oli
|