This subroutine calculates the gradient of the energy, gb, but calls the subroutine GBAH to do so. In the serial version, the majority of the subroutine ENERGY was dedicated to calculating the gb array. It was much easier to create a subroutine for this calculation since the next subroutine, DEMBX_MPI, has the same functional form for creating its Ah array.
The differences here appear when calculating utot. One needs to get the per-processor value dutot and then make a global sum with a call to MPI_ALLREDUCE. The last two lines in the subroutine are simple assignments:
utot = utot + C
gb = gb + b
The first line does an update on a single variable, utot. But the second line does an update of the entire array gb with array b using FORTRAN90 syntax. Keep in mind, each processor is just updating its own local copy of gb with its local copy of b.
In THERMAL3D_MPI.f, there is additional complexity for this subroutine. After it updates gb (like above), it must calculate its value for the "constant" macrostrain energy term, C, using the zcon and u2 arrays. This value C is added to utot (also above), but to get the final value of utot, one must add the constant Y from FEMAT_MPI.