When calculating the gb or Ah arrays in the serial program, one can see that they are essentially the same calculation with a 1:1 correspondence between the arrays u and h as well as gb and Ah. Therefore, it is logical to create this subroutine and give input parameters which determine if one is calculating gb or Ah. Regardless if one wants gb or Ah, the local array that is calculated in GBAH is the array om; the place from which GBAH is called and its parameters determine if gb or Ah is going to be calculated.
In this subroutine, elements of the array om (output matrix) are being calculated on a per node basis. In addition to that, FORTRAN90 syntax to sum over all values of an array index to calculate the terms removes the extra loop (n = 1, 3) that was found in the serial version. An example calculation from ELAS3D_MPI.f and ELAS3D.f should illustrate the point. Note that the last index of om, j, would not appear in ELECFEM3D_MPI.f (cf. Table 1).
The u(ib(m, 1), n) and u(ib(m, 2), n) terms in the comments refer to the like terms found in the serial version, ie.
om exhibits the same behavior as Ah and gb, in other words, it must also have top and bottom ghost layers as well. Therefore, before GBAH returns om back to the calling routine, it has already created the required ghost layers by judicious use of the subroutines t2b_dp and b2t_dp.
But the most intriguing part of the calculation involves the determination for the necessary elements of uh and vox in the above calculation. Note that uh is merely a parameter; uh equals u when called from subroutine ENERGY_MPI but equals h when called from DEMBX_MPI. This parameter is named uh to remind the user of its duality.
Before the code was written, it was necessary to determine what the (i, j, k) indices for the uh and vox arrays for each of the 27 terms will be, keeping in mind ib and Table 3 in the serial code manual (NISTIR 6269). After one has done this, their (i, j, k) indices can be deduced. However, one must also keep in mind the periodic nature of uh and vox in the x and y directions; periodicity in the z direction is assured due to the ghost layers. Therefore when looping over im and jm (ie. x and y directions), special precautions (in the form of if-statements) are taken when (im + 1) > nx or (im − 1) <, when the value of that particular i-type coordinate is assigned to 1 or nx, respectively, with similar arrangements for jm.
For a given set of (im, jm, km) loop variables, if the serial code says to invoke ib(m, 1), this corresponds to the listing in Table 3 in the serial code manual (NISTIR 6269) for ib(m, 1), which says the required neighbor has a (Δi,Δj,Δk) setting of (0,1,0). This implies that the proper element to use in this calculation would be (im, jm + 1, km). There is a change in variable for jm + 1 which is now assigned to ifya. This specialized nomenclature is summarized here:
Therefore, when one encounters a ifxa, ifxb, ifya or ifyb term, it can immediately be identified as a term that needs to be fixed or adjusted according to the above specifications. This occurs not only for uh, but also for vox(im, jm, km) as it is a parameter for dk.