This subroutine calculates the gradient relaxation process and the h and Ah arrays, but calls the subroutine GBAH to generate them. It behaves similarly to the serial version, except in cases when FORTRAN90 array syntax was implemented for lines like:
| u=u-lambda*h |
| gb=gb-lambda*Ah |
In the FORTRAN77 serial version, these lines looked like:
| do 540 m3=1,3 | |
| do 540 m=1,ns | |
| u(m,m3)=u(m,m3)-lambda*h(m,m3) | |
| gb(m,m3)=gb(m,m3)-lambda*Ah(m,m3) | |
| 540 | continue |