39  module subroutine opr_cpu_opgrad(ux, uy, uz, u, coef, e_start, e_end)
 
   40    type(coef_t), 
intent(in) :: coef
 
   41    integer, 
intent(in) :: e_start, e_end
 
   42    real(kind=rp), 
dimension(coef%Xh%lxyz, e_end - e_start + 1), &
 
   44    real(kind=rp), 
dimension(coef%Xh%lxyz, e_end - e_start + 1), &
 
   46    real(kind=rp), 
dimension(coef%Xh%lxyz, e_end - e_start + 1), &
 
   48    real(kind=rp), 
dimension(coef%Xh%lxyz, e_end - e_start + 1), &
 
   52    e_len = e_end-e_start+1
 
   54    if (e_len .eq. 1) 
then 
   55       call opr_cpu_opgrad_single(ux, uy, uz, u, coef, e_start)
 
   57       call opr_cpu_opgrad_many(ux, uy, uz, u, coef, e_start, e_len)
 
   59  end subroutine opr_cpu_opgrad
 
   61  subroutine opr_cpu_opgrad_many(ux, uy, uz, u, coef, e_start, e_len)
 
   62    type(coef_t), 
intent(in) :: coef
 
   63    integer, 
intent(in) :: e_start, e_len
 
   64    real(kind=rp), 
dimension(coef%Xh%lxyz, e_len), 
intent(inout) :: ux
 
   65    real(kind=rp), 
dimension(coef%Xh%lxyz, e_len), 
intent(inout) :: uy
 
   66    real(kind=rp), 
dimension(coef%Xh%lxyz, e_len), 
intent(inout) :: uz
 
   67    real(kind=rp), 
dimension(coef%Xh%lxyz, e_len), 
intent(in) :: u
 
   69    associate(xh => coef%Xh, msh => coef%msh, &
 
   70         drdx => coef%drdx, drdy => coef%drdy, drdz => coef%drdz, &
 
   71         dsdx => coef%dsdx, dsdy => coef%dsdy, dsdz => coef%dsdz, &
 
   72         dtdx => coef%dtdx, dtdy => coef%dtdy, dtdz => coef%dtdz)
 
   76         call cpu_opgrad_lx18(ux, uy, uz, u, &
 
   77              xh%dx, xh%dy, xh%dz, &
 
   78              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
   79              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
   80              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
   83         call cpu_opgrad_lx17(ux, uy, uz, u, &
 
   84              xh%dx, xh%dy, xh%dz, &
 
   85              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
   86              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
   87              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
   90         call cpu_opgrad_lx16(ux, uy, uz, u, &
 
   91              xh%dx, xh%dy, xh%dz, &
 
   92              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
   93              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
   94              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
   97         call cpu_opgrad_lx15(ux, uy, uz, u, &
 
   98              xh%dx, xh%dy, xh%dz, &
 
   99              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  100              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  101              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  104         call cpu_opgrad_lx14(ux, uy, uz, u, &
 
  105              xh%dx, xh%dy, xh%dz, &
 
  106              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  107              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  108              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  111         call cpu_opgrad_lx13(ux, uy, uz, u, &
 
  112              xh%dx, xh%dy, xh%dz, &
 
  113              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  114              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  115              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  118         call cpu_opgrad_lx12(ux, uy, uz, u, &
 
  119              xh%dx, xh%dy, xh%dz, &
 
  120              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  121              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  122              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  125         call cpu_opgrad_lx11(ux, uy, uz, u, &
 
  126              xh%dx, xh%dy, xh%dz, &
 
  127              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  128              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  129              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  132         call cpu_opgrad_lx10(ux, uy, uz, u, &
 
  133              xh%dx, xh%dy, xh%dz, &
 
  134              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  135              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  136              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  140         call cpu_opgrad_lx9(ux, uy, uz, u, &
 
  141              xh%dx, xh%dy, xh%dz, &
 
  142              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  143              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  144              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  147         call cpu_opgrad_lx8(ux, uy, uz, u, &
 
  148              xh%dx, xh%dy, xh%dz, &
 
  149              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  150              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  151              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  154         call cpu_opgrad_lx7(ux, uy, uz, u, &
 
  155              xh%dx, xh%dy, xh%dz, &
 
  156              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  157              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  158              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  161         call cpu_opgrad_lx6(ux, uy, uz, u, &
 
  162              xh%dx, xh%dy, xh%dz, &
 
  163              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  164              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  165              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  168         call cpu_opgrad_lx5(ux, uy, uz, u, &
 
  169              xh%dx, xh%dy, xh%dz, &
 
  170              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  171              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  172              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  175         call cpu_opgrad_lx4(ux, uy, uz, u, &
 
  176              xh%dx, xh%dy, xh%dz, &
 
  177              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  178              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  179              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  182         call cpu_opgrad_lx3(ux, uy, uz, u, &
 
  183              xh%dx, xh%dy, xh%dz, &
 
  184              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  185              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  186              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  189         call cpu_opgrad_lx2(ux, uy, uz, u, &
 
  190              xh%dx, xh%dy, xh%dz, &
 
  191              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  192              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  193              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  196         call cpu_opgrad_lx(ux, uy, uz, u, &
 
  197              xh%dx, xh%dy, xh%dz, &
 
  198              drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
 
  199              drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
 
  200              drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
 
  205  end subroutine opr_cpu_opgrad_many
 
  207  subroutine cpu_opgrad_lx(ux, uy, uz, u, dx, dy, dz, &
 
  208       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n, lx)
 
  209    integer, 
intent(in) :: n, lx
 
  210    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
  211    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
  212    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
  213    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
  214    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
  215    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
  216    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
  217    real(kind=rp) :: ur(lx, lx, lx)
 
  218    real(kind=rp) :: us(lx, lx, lx)
 
  219    real(kind=rp) :: ut(lx, lx, lx)
 
  221    integer :: e, i, j, k, l
 
  228                tmp = tmp + dx(i,k) * u(k,j,1,e)
 
  239                   tmp = tmp + dy(j,l) * u(i,l,k,e)
 
  250                tmp = tmp + dz(k,l) * u(i,1,l,e)
 
  256       do i = 1, lx * lx * lx
 
  257          ux(i,1,1,e) = w3(i,1,1) &
 
  258                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
  259                       + dsdx(i,1,1,e) * us(i,1,1) &
 
  260                       + dtdx(i,1,1,e) * ut(i,1,1) )
 
  261          uy(i,1,1,e) = w3(i,1,1) &
 
  262                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
  263                        + drdy(i,1,1,e) * ur(i,1,1) &
 
  264                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
  265          uz(i,1,1,e) = w3(i,1,1) &
 
  266                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
  267                        + drdz(i,1,1,e) * ur(i,1,1) &
 
  268                        + dsdz(i,1,1,e) * us(i,1,1) )
 
  271  end subroutine cpu_opgrad_lx
 
  273  subroutine cpu_opgrad_lx18(ux, uy, uz, u, dx, dy, dz, &
 
  274       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
  275    integer, 
parameter :: lx = 18
 
  276    integer, 
intent(in) :: n
 
  277    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
  278    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
  279    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
  280    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
  281    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
  282    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
  283    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
  284    real(kind=rp) :: ur(lx, lx, lx)
 
  285    real(kind=rp) :: us(lx, lx, lx)
 
  286    real(kind=rp) :: ut(lx, lx, lx)
 
  287    integer :: e, i, j, k
 
  292             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
  293                        + dx(i,2) * u(2,j,1,e) &
 
  294                        + dx(i,3) * u(3,j,1,e) &
 
  295                        + dx(i,4) * u(4,j,1,e) &
 
  296                        + dx(i,5) * u(5,j,1,e) &
 
  297                        + dx(i,6) * u(6,j,1,e) &
 
  298                        + dx(i,7) * u(7,j,1,e) &
 
  299                        + dx(i,8) * u(8,j,1,e) &
 
  300                        + dx(i,9) * u(9,j,1,e) &
 
  301                        + dx(i,10) * u(10,j,1,e) &
 
  302                        + dx(i,11) * u(11,j,1,e) &
 
  303                        + dx(i,12) * u(12,j,1,e) &
 
  304                        + dx(i,13) * u(13,j,1,e) &
 
  305                        + dx(i,14) * u(14,j,1,e) &
 
  306                        + dx(i,15) * u(15,j,1,e) &
 
  307                        + dx(i,16) * u(16,j,1,e) &
 
  308                        + dx(i,17) * u(17,j,1,e) &
 
  309                        + dx(i,18) * u(18,j,1,e)
 
  316                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
  317                          + dy(j,2) * u(i,2,k,e) &
 
  318                          + dy(j,3) * u(i,3,k,e) &
 
  319                          + dy(j,4) * u(i,4,k,e) &
 
  320                          + dy(j,5) * u(i,5,k,e) &
 
  321                          + dy(j,6) * u(i,6,k,e) &
 
  322                          + dy(j,7) * u(i,7,k,e) &
 
  323                          + dy(j,8) * u(i,8,k,e) &
 
  324                          + dy(j,9) * u(i,9,k,e) &
 
  325                          + dy(j,10) * u(i,10,k,e) &
 
  326                          + dy(j,11) * u(i,11,k,e) &
 
  327                          + dy(j,12) * u(i,12,k,e) &
 
  328                          + dy(j,13) * u(i,13,k,e) &
 
  329                          + dy(j,14) * u(i,14,k,e) &
 
  330                          + dy(j,15) * u(i,15,k,e) &
 
  331                          + dy(j,16) * u(i,16,k,e) &
 
  332                          + dy(j,17) * u(i,17,k,e) &
 
  333                          + dy(j,18) * u(i,18,k,e)
 
  340             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
  341                       + dz(k,2) * u(i,1,2,e) &
 
  342                       + dz(k,3) * u(i,1,3,e) &
 
  343                       + dz(k,4) * u(i,1,4,e) &
 
  344                       + dz(k,5) * u(i,1,5,e) &
 
  345                       + dz(k,6) * u(i,1,6,e) &
 
  346                       + dz(k,7) * u(i,1,7,e) &
 
  347                       + dz(k,8) * u(i,1,8,e) &
 
  348                       + dz(k,9) * u(i,1,9,e) &
 
  349                       + dz(k,10) * u(i,1,10,e) &
 
  350                       + dz(k,11) * u(i,1,11,e) &
 
  351                       + dz(k,12) * u(i,1,12,e) &
 
  352                       + dz(k,13) * u(i,1,13,e) &
 
  353                       + dz(k,14) * u(i,1,14,e) &
 
  354                       + dz(k,15) * u(i,1,15,e) &
 
  355                       + dz(k,16) * u(i,1,16,e) &
 
  356                       + dz(k,17) * u(i,1,17,e) &
 
  357                       + dz(k,18) * u(i,1,18,e)
 
  361       do i = 1, lx * lx * lx
 
  362          ux(i,1,1,e) = w3(i,1,1) &
 
  363                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
  364                       + dsdx(i,1,1,e) * us(i,1,1) &
 
  365                       + dtdx(i,1,1,e) * ut(i,1,1) )
 
  366          uy(i,1,1,e) = w3(i,1,1) &
 
  367                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
  368                        + drdy(i,1,1,e) * ur(i,1,1) &
 
  369                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
  370          uz(i,1,1,e) = w3(i,1,1) &
 
  371                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
  372                        + drdz(i,1,1,e) * ur(i,1,1) &
 
  373                        + dsdz(i,1,1,e) * us(i,1,1) )
 
  376  end subroutine cpu_opgrad_lx18
 
  378  subroutine cpu_opgrad_lx17(ux, uy, uz, u, dx, dy, dz, &
 
  379       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
  380    integer, 
parameter :: lx = 17
 
  381    integer, 
intent(in) :: n
 
  382    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
  383    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
  384    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
  385    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
  386    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
  387    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
  388    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
  389    real(kind=rp) :: ur(lx, lx, lx)
 
  390    real(kind=rp) :: us(lx, lx, lx)
 
  391    real(kind=rp) :: ut(lx, lx, lx)
 
  392    integer :: e, i, j, k
 
  397             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
  398                        + dx(i,2) * u(2,j,1,e) &
 
  399                        + dx(i,3) * u(3,j,1,e) &
 
  400                        + dx(i,4) * u(4,j,1,e) &
 
  401                        + dx(i,5) * u(5,j,1,e) &
 
  402                        + dx(i,6) * u(6,j,1,e) &
 
  403                        + dx(i,7) * u(7,j,1,e) &
 
  404                        + dx(i,8) * u(8,j,1,e) &
 
  405                        + dx(i,9) * u(9,j,1,e) &
 
  406                        + dx(i,10) * u(10,j,1,e) &
 
  407                        + dx(i,11) * u(11,j,1,e) &
 
  408                        + dx(i,12) * u(12,j,1,e) &
 
  409                        + dx(i,13) * u(13,j,1,e) &
 
  410                        + dx(i,14) * u(14,j,1,e) &
 
  411                        + dx(i,15) * u(15,j,1,e) &
 
  412                        + dx(i,16) * u(16,j,1,e) &
 
  413                        + dx(i,17) * u(17,j,1,e)
 
  420                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
  421                          + dy(j,2) * u(i,2,k,e) &
 
  422                          + dy(j,3) * u(i,3,k,e) &
 
  423                          + dy(j,4) * u(i,4,k,e) &
 
  424                          + dy(j,5) * u(i,5,k,e) &
 
  425                          + dy(j,6) * u(i,6,k,e) &
 
  426                          + dy(j,7) * u(i,7,k,e) &
 
  427                          + dy(j,8) * u(i,8,k,e) &
 
  428                          + dy(j,9) * u(i,9,k,e) &
 
  429                          + dy(j,10) * u(i,10,k,e) &
 
  430                          + dy(j,11) * u(i,11,k,e) &
 
  431                          + dy(j,12) * u(i,12,k,e) &
 
  432                          + dy(j,13) * u(i,13,k,e) &
 
  433                          + dy(j,14) * u(i,14,k,e) &
 
  434                          + dy(j,15) * u(i,15,k,e) &
 
  435                          + dy(j,16) * u(i,16,k,e) &
 
  436                          + dy(j,17) * u(i,17,k,e)
 
  443             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
  444                       + dz(k,2) * u(i,1,2,e) &
 
  445                       + dz(k,3) * u(i,1,3,e) &
 
  446                       + dz(k,4) * u(i,1,4,e) &
 
  447                       + dz(k,5) * u(i,1,5,e) &
 
  448                       + dz(k,6) * u(i,1,6,e) &
 
  449                       + dz(k,7) * u(i,1,7,e) &
 
  450                       + dz(k,8) * u(i,1,8,e) &
 
  451                       + dz(k,9) * u(i,1,9,e) &
 
  452                       + dz(k,10) * u(i,1,10,e) &
 
  453                       + dz(k,11) * u(i,1,11,e) &
 
  454                       + dz(k,12) * u(i,1,12,e) &
 
  455                       + dz(k,13) * u(i,1,13,e) &
 
  456                       + dz(k,14) * u(i,1,14,e) &
 
  457                       + dz(k,15) * u(i,1,15,e) &
 
  458                       + dz(k,16) * u(i,1,16,e) &
 
  459                       + dz(k,17) * u(i,1,17,e)
 
  463       do i = 1, lx * lx * lx
 
  464          ux(i,1,1,e) = w3(i,1,1) &
 
  465                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
  466                       + dsdx(i,1,1,e) * us(i,1,1) &
 
  467                       + dtdx(i,1,1,e) * ut(i,1,1) )
 
  468          uy(i,1,1,e) = w3(i,1,1) &
 
  469                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
  470                        + drdy(i,1,1,e) * ur(i,1,1) &
 
  471                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
  472          uz(i,1,1,e) = w3(i,1,1) &
 
  473                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
  474                        + drdz(i,1,1,e) * ur(i,1,1) &
 
  475                        + dsdz(i,1,1,e) * us(i,1,1) )
 
  478  end subroutine cpu_opgrad_lx17
 
  480  subroutine cpu_opgrad_lx16(ux, uy, uz, u, dx, dy, dz, &
 
  481       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
  482    integer, 
parameter :: lx = 16
 
  483    integer, 
intent(in) :: n
 
  484    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
  485    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
  486    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
  487    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
  488    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
  489    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
  490    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
  491    real(kind=rp) :: ur(lx, lx, lx)
 
  492    real(kind=rp) :: us(lx, lx, lx)
 
  493    real(kind=rp) :: ut(lx, lx, lx)
 
  494    integer :: e, i, j, k
 
  499             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
  500                        + dx(i,2) * u(2,j,1,e) &
 
  501                        + dx(i,3) * u(3,j,1,e) &
 
  502                        + dx(i,4) * u(4,j,1,e) &
 
  503                        + dx(i,5) * u(5,j,1,e) &
 
  504                        + dx(i,6) * u(6,j,1,e) &
 
  505                        + dx(i,7) * u(7,j,1,e) &
 
  506                        + dx(i,8) * u(8,j,1,e) &
 
  507                        + dx(i,9) * u(9,j,1,e) &
 
  508                        + dx(i,10) * u(10,j,1,e) &
 
  509                        + dx(i,11) * u(11,j,1,e) &
 
  510                        + dx(i,12) * u(12,j,1,e) &
 
  511                        + dx(i,13) * u(13,j,1,e) &
 
  512                        + dx(i,14) * u(14,j,1,e) &
 
  513                        + dx(i,15) * u(15,j,1,e) &
 
  514                        + dx(i,16) * u(16,j,1,e)
 
  521                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
  522                          + dy(j,2) * u(i,2,k,e) &
 
  523                          + dy(j,3) * u(i,3,k,e) &
 
  524                          + dy(j,4) * u(i,4,k,e) &
 
  525                          + dy(j,5) * u(i,5,k,e) &
 
  526                          + dy(j,6) * u(i,6,k,e) &
 
  527                          + dy(j,7) * u(i,7,k,e) &
 
  528                          + dy(j,8) * u(i,8,k,e) &
 
  529                          + dy(j,9) * u(i,9,k,e) &
 
  530                          + dy(j,10) * u(i,10,k,e) &
 
  531                          + dy(j,11) * u(i,11,k,e) &
 
  532                          + dy(j,12) * u(i,12,k,e) &
 
  533                          + dy(j,13) * u(i,13,k,e) &
 
  534                          + dy(j,14) * u(i,14,k,e) &
 
  535                          + dy(j,15) * u(i,15,k,e) &
 
  536                          + dy(j,16) * u(i,16,k,e)
 
  543             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
  544                       + dz(k,2) * u(i,1,2,e) &
 
  545                       + dz(k,3) * u(i,1,3,e) &
 
  546                       + dz(k,4) * u(i,1,4,e) &
 
  547                       + dz(k,5) * u(i,1,5,e) &
 
  548                       + dz(k,6) * u(i,1,6,e) &
 
  549                       + dz(k,7) * u(i,1,7,e) &
 
  550                       + dz(k,8) * u(i,1,8,e) &
 
  551                       + dz(k,9) * u(i,1,9,e) &
 
  552                       + dz(k,10) * u(i,1,10,e) &
 
  553                       + dz(k,11) * u(i,1,11,e) &
 
  554                       + dz(k,12) * u(i,1,12,e) &
 
  555                       + dz(k,13) * u(i,1,13,e) &
 
  556                       + dz(k,14) * u(i,1,14,e) &
 
  557                       + dz(k,15) * u(i,1,15,e) &
 
  558                       + dz(k,16) * u(i,1,16,e)
 
  562       do i = 1, lx * lx * lx
 
  563          ux(i,1,1,e) = w3(i,1,1) &
 
  564                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
  565                       + dsdx(i,1,1,e) * us(i,1,1) &
 
  566                       + dtdx(i,1,1,e) * ut(i,1,1) )
 
  567          uy(i,1,1,e) = w3(i,1,1) &
 
  568                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
  569                        + drdy(i,1,1,e) * ur(i,1,1) &
 
  570                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
  571          uz(i,1,1,e) = w3(i,1,1) &
 
  572                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
  573                        + drdz(i,1,1,e) * ur(i,1,1) &
 
  574                        + dsdz(i,1,1,e) * us(i,1,1) )
 
  577  end subroutine cpu_opgrad_lx16
 
  579  subroutine cpu_opgrad_lx15(ux, uy, uz, u, dx, dy, dz, &
 
  580       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
  581    integer, 
parameter :: lx = 15
 
  582    integer, 
intent(in) :: n
 
  583    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
  584    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
  585    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
  586    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
  587    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
  588    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
  589    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
  590    real(kind=rp) :: ur(lx, lx, lx)
 
  591    real(kind=rp) :: us(lx, lx, lx)
 
  592    real(kind=rp) :: ut(lx, lx, lx)
 
  593    integer :: e, i, j, k
 
  598             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
  599                        + dx(i,2) * u(2,j,1,e) &
 
  600                        + dx(i,3) * u(3,j,1,e) &
 
  601                        + dx(i,4) * u(4,j,1,e) &
 
  602                        + dx(i,5) * u(5,j,1,e) &
 
  603                        + dx(i,6) * u(6,j,1,e) &
 
  604                        + dx(i,7) * u(7,j,1,e) &
 
  605                        + dx(i,8) * u(8,j,1,e) &
 
  606                        + dx(i,9) * u(9,j,1,e) &
 
  607                        + dx(i,10) * u(10,j,1,e) &
 
  608                        + dx(i,11) * u(11,j,1,e) &
 
  609                        + dx(i,12) * u(12,j,1,e) &
 
  610                        + dx(i,13) * u(13,j,1,e) &
 
  611                        + dx(i,14) * u(14,j,1,e) &
 
  612                        + dx(i,15) * u(15,j,1,e)
 
  619                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
  620                          + dy(j,2) * u(i,2,k,e) &
 
  621                          + dy(j,3) * u(i,3,k,e) &
 
  622                          + dy(j,4) * u(i,4,k,e) &
 
  623                          + dy(j,5) * u(i,5,k,e) &
 
  624                          + dy(j,6) * u(i,6,k,e) &
 
  625                          + dy(j,7) * u(i,7,k,e) &
 
  626                          + dy(j,8) * u(i,8,k,e) &
 
  627                          + dy(j,9) * u(i,9,k,e) &
 
  628                          + dy(j,10) * u(i,10,k,e) &
 
  629                          + dy(j,11) * u(i,11,k,e) &
 
  630                          + dy(j,12) * u(i,12,k,e) &
 
  631                          + dy(j,13) * u(i,13,k,e) &
 
  632                          + dy(j,14) * u(i,14,k,e) &
 
  633                          + dy(j,15) * u(i,15,k,e)
 
  640             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
  641                       + dz(k,2) * u(i,1,2,e) &
 
  642                       + dz(k,3) * u(i,1,3,e) &
 
  643                       + dz(k,4) * u(i,1,4,e) &
 
  644                       + dz(k,5) * u(i,1,5,e) &
 
  645                       + dz(k,6) * u(i,1,6,e) &
 
  646                       + dz(k,7) * u(i,1,7,e) &
 
  647                       + dz(k,8) * u(i,1,8,e) &
 
  648                       + dz(k,9) * u(i,1,9,e) &
 
  649                       + dz(k,10) * u(i,1,10,e) &
 
  650                       + dz(k,11) * u(i,1,11,e) &
 
  651                       + dz(k,12) * u(i,1,12,e) &
 
  652                       + dz(k,13) * u(i,1,13,e) &
 
  653                       + dz(k,14) * u(i,1,14,e) &
 
  654                       + dz(k,15) * u(i,1,15,e)
 
  658       do i = 1, lx * lx * lx
 
  659          ux(i,1,1,e) = w3(i,1,1) &
 
  660                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
  661                       + dsdx(i,1,1,e) * us(i,1,1) &
 
  662                       + dtdx(i,1,1,e) * ut(i,1,1) )
 
  663          uy(i,1,1,e) = w3(i,1,1) &
 
  664                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
  665                        + drdy(i,1,1,e) * ur(i,1,1) &
 
  666                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
  667          uz(i,1,1,e) = w3(i,1,1) &
 
  668                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
  669                        + drdz(i,1,1,e) * ur(i,1,1) &
 
  670                        + dsdz(i,1,1,e) * us(i,1,1) )
 
  673  end subroutine cpu_opgrad_lx15
 
  675  subroutine cpu_opgrad_lx14(ux, uy, uz, u, dx, dy, dz, &
 
  676       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
  677    integer, 
parameter :: lx = 14
 
  678    integer, 
intent(in) :: n
 
  679    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
  680    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
  681    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
  682    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
  683    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
  684    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
  685    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
  686    real(kind=rp) :: ur(lx, lx, lx)
 
  687    real(kind=rp) :: us(lx, lx, lx)
 
  688    real(kind=rp) :: ut(lx, lx, lx)
 
  689    integer :: e, i, j, k
 
  694             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
  695                        + dx(i,2) * u(2,j,1,e) &
 
  696                        + dx(i,3) * u(3,j,1,e) &
 
  697                        + dx(i,4) * u(4,j,1,e) &
 
  698                        + dx(i,5) * u(5,j,1,e) &
 
  699                        + dx(i,6) * u(6,j,1,e) &
 
  700                        + dx(i,7) * u(7,j,1,e) &
 
  701                        + dx(i,8) * u(8,j,1,e) &
 
  702                        + dx(i,9) * u(9,j,1,e) &
 
  703                        + dx(i,10) * u(10,j,1,e) &
 
  704                        + dx(i,11) * u(11,j,1,e) &
 
  705                        + dx(i,12) * u(12,j,1,e) &
 
  706                        + dx(i,13) * u(13,j,1,e) &
 
  707                        + dx(i,14) * u(14,j,1,e)
 
  714                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
  715                          + dy(j,2) * u(i,2,k,e) &
 
  716                          + dy(j,3) * u(i,3,k,e) &
 
  717                          + dy(j,4) * u(i,4,k,e) &
 
  718                          + dy(j,5) * u(i,5,k,e) &
 
  719                          + dy(j,6) * u(i,6,k,e) &
 
  720                          + dy(j,7) * u(i,7,k,e) &
 
  721                          + dy(j,8) * u(i,8,k,e) &
 
  722                          + dy(j,9) * u(i,9,k,e) &
 
  723                          + dy(j,10) * u(i,10,k,e) &
 
  724                          + dy(j,11) * u(i,11,k,e) &
 
  725                          + dy(j,12) * u(i,12,k,e) &
 
  726                          + dy(j,13) * u(i,13,k,e) &
 
  727                          + dy(j,14) * u(i,14,k,e)
 
  734             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
  735                       + dz(k,2) * u(i,1,2,e) &
 
  736                       + dz(k,3) * u(i,1,3,e) &
 
  737                       + dz(k,4) * u(i,1,4,e) &
 
  738                       + dz(k,5) * u(i,1,5,e) &
 
  739                       + dz(k,6) * u(i,1,6,e) &
 
  740                       + dz(k,7) * u(i,1,7,e) &
 
  741                       + dz(k,8) * u(i,1,8,e) &
 
  742                       + dz(k,9) * u(i,1,9,e) &
 
  743                       + dz(k,10) * u(i,1,10,e) &
 
  744                       + dz(k,11) * u(i,1,11,e) &
 
  745                       + dz(k,12) * u(i,1,12,e) &
 
  746                       + dz(k,13) * u(i,1,13,e) &
 
  747                       + dz(k,14) * u(i,1,14,e)
 
  751       do i = 1, lx * lx * lx
 
  752          ux(i,1,1,e) = w3(i,1,1) &
 
  753                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
  754                       + dsdx(i,1,1,e) * us(i,1,1) &
 
  755                       + dtdx(i,1,1,e) * ut(i,1,1) )
 
  756          uy(i,1,1,e) = w3(i,1,1) &
 
  757                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
  758                        + drdy(i,1,1,e) * ur(i,1,1) &
 
  759                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
  760          uz(i,1,1,e) = w3(i,1,1) &
 
  761                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
  762                        + drdz(i,1,1,e) * ur(i,1,1) &
 
  763                        + dsdz(i,1,1,e) * us(i,1,1) )
 
  766  end subroutine cpu_opgrad_lx14
 
  768  subroutine cpu_opgrad_lx13(ux, uy, uz, u, dx, dy, dz, &
 
  769       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
  770    integer, 
parameter :: lx = 13
 
  771    integer, 
intent(in) :: n
 
  772    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
  773    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
  774    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
  775    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
  776    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
  777    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
  778    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
  779    real(kind=rp) :: ur(lx, lx, lx)
 
  780    real(kind=rp) :: us(lx, lx, lx)
 
  781    real(kind=rp) :: ut(lx, lx, lx)
 
  782    integer :: e, i, j, k
 
  787             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
  788                        + dx(i,2) * u(2,j,1,e) &
 
  789                        + dx(i,3) * u(3,j,1,e) &
 
  790                        + dx(i,4) * u(4,j,1,e) &
 
  791                        + dx(i,5) * u(5,j,1,e) &
 
  792                        + dx(i,6) * u(6,j,1,e) &
 
  793                        + dx(i,7) * u(7,j,1,e) &
 
  794                        + dx(i,8) * u(8,j,1,e) &
 
  795                        + dx(i,9) * u(9,j,1,e) &
 
  796                        + dx(i,10) * u(10,j,1,e) &
 
  797                        + dx(i,11) * u(11,j,1,e) &
 
  798                        + dx(i,12) * u(12,j,1,e) &
 
  799                        + dx(i,13) * u(13,j,1,e)
 
  806                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
  807                          + dy(j,2) * u(i,2,k,e) &
 
  808                          + dy(j,3) * u(i,3,k,e) &
 
  809                          + dy(j,4) * u(i,4,k,e) &
 
  810                          + dy(j,5) * u(i,5,k,e) &
 
  811                          + dy(j,6) * u(i,6,k,e) &
 
  812                          + dy(j,7) * u(i,7,k,e) &
 
  813                          + dy(j,8) * u(i,8,k,e) &
 
  814                          + dy(j,9) * u(i,9,k,e) &
 
  815                          + dy(j,10) * u(i,10,k,e) &
 
  816                          + dy(j,11) * u(i,11,k,e) &
 
  817                          + dy(j,12) * u(i,12,k,e) &
 
  818                          + dy(j,13) * u(i,13,k,e)
 
  825             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
  826                       + dz(k,2) * u(i,1,2,e) &
 
  827                       + dz(k,3) * u(i,1,3,e) &
 
  828                       + dz(k,4) * u(i,1,4,e) &
 
  829                       + dz(k,5) * u(i,1,5,e) &
 
  830                       + dz(k,6) * u(i,1,6,e) &
 
  831                       + dz(k,7) * u(i,1,7,e) &
 
  832                       + dz(k,8) * u(i,1,8,e) &
 
  833                       + dz(k,9) * u(i,1,9,e) &
 
  834                       + dz(k,10) * u(i,1,10,e) &
 
  835                       + dz(k,11) * u(i,1,11,e) &
 
  836                       + dz(k,12) * u(i,1,12,e) &
 
  837                       + dz(k,13) * u(i,1,13,e)
 
  841       do i = 1, lx * lx * lx
 
  842          ux(i,1,1,e) = w3(i,1,1) &
 
  843                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
  844                       + dsdx(i,1,1,e) * us(i,1,1) &
 
  845                       + dtdx(i,1,1,e) * ut(i,1,1) )
 
  846          uy(i,1,1,e) = w3(i,1,1) &
 
  847                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
  848                        + drdy(i,1,1,e) * ur(i,1,1) &
 
  849                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
  850          uz(i,1,1,e) = w3(i,1,1) &
 
  851                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
  852                        + drdz(i,1,1,e) * ur(i,1,1) &
 
  853                        + dsdz(i,1,1,e) * us(i,1,1) )
 
  856  end subroutine cpu_opgrad_lx13
 
  858  subroutine cpu_opgrad_lx12(ux, uy, uz, u, dx, dy, dz, &
 
  859       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
  860    integer, 
parameter :: lx = 12
 
  861    integer, 
intent(in) :: n
 
  862    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
  863    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
  864    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
  865    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
  866    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
  867    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
  868    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
  869    real(kind=rp) :: ur(lx, lx, lx)
 
  870    real(kind=rp) :: us(lx, lx, lx)
 
  871    real(kind=rp) :: ut(lx, lx, lx)
 
  872    integer :: e, i, j, k
 
  877             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
  878                        + dx(i,2) * u(2,j,1,e) &
 
  879                        + dx(i,3) * u(3,j,1,e) &
 
  880                        + dx(i,4) * u(4,j,1,e) &
 
  881                        + dx(i,5) * u(5,j,1,e) &
 
  882                        + dx(i,6) * u(6,j,1,e) &
 
  883                        + dx(i,7) * u(7,j,1,e) &
 
  884                        + dx(i,8) * u(8,j,1,e) &
 
  885                        + dx(i,9) * u(9,j,1,e) &
 
  886                        + dx(i,10) * u(10,j,1,e) &
 
  887                        + dx(i,11) * u(11,j,1,e) &
 
  888                        + dx(i,12) * u(12,j,1,e)
 
  895                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
  896                          + dy(j,2) * u(i,2,k,e) &
 
  897                          + dy(j,3) * u(i,3,k,e) &
 
  898                          + dy(j,4) * u(i,4,k,e) &
 
  899                          + dy(j,5) * u(i,5,k,e) &
 
  900                          + dy(j,6) * u(i,6,k,e) &
 
  901                          + dy(j,7) * u(i,7,k,e) &
 
  902                          + dy(j,8) * u(i,8,k,e) &
 
  903                          + dy(j,9) * u(i,9,k,e) &
 
  904                          + dy(j,10) * u(i,10,k,e) &
 
  905                          + dy(j,11) * u(i,11,k,e) &
 
  906                          + dy(j,12) * u(i,12,k,e)
 
  913             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
  914                       + dz(k,2) * u(i,1,2,e) &
 
  915                       + dz(k,3) * u(i,1,3,e) &
 
  916                       + dz(k,4) * u(i,1,4,e) &
 
  917                       + dz(k,5) * u(i,1,5,e) &
 
  918                       + dz(k,6) * u(i,1,6,e) &
 
  919                       + dz(k,7) * u(i,1,7,e) &
 
  920                       + dz(k,8) * u(i,1,8,e) &
 
  921                       + dz(k,9) * u(i,1,9,e) &
 
  922                       + dz(k,10) * u(i,1,10,e) &
 
  923                       + dz(k,11) * u(i,1,11,e) &
 
  924                       + dz(k,12) * u(i,1,12,e)
 
  928       do i = 1, lx * lx * lx
 
  929          ux(i,1,1,e) = w3(i,1,1) &
 
  930                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
  931                       + dsdx(i,1,1,e) * us(i,1,1) &
 
  932                       + dtdx(i,1,1,e) * ut(i,1,1) )
 
  933          uy(i,1,1,e) = w3(i,1,1) &
 
  934                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
  935                        + drdy(i,1,1,e) * ur(i,1,1) &
 
  936                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
  937          uz(i,1,1,e) = w3(i,1,1) &
 
  938                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
  939                        + drdz(i,1,1,e) * ur(i,1,1) &
 
  940                        + dsdz(i,1,1,e) * us(i,1,1) )
 
  943  end subroutine cpu_opgrad_lx12
 
  945  subroutine cpu_opgrad_lx11(ux, uy, uz, u, dx, dy, dz, &
 
  946       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
  947    integer, 
parameter :: lx = 11
 
  948    integer, 
intent(in) :: n
 
  949    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
  950    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
  951    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
  952    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
  953    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
  954    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
  955    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
  956    real(kind=rp) :: ur(lx, lx, lx)
 
  957    real(kind=rp) :: us(lx, lx, lx)
 
  958    real(kind=rp) :: ut(lx, lx, lx)
 
  959    integer :: e, i, j, k
 
  964             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
  965                        + dx(i,2) * u(2,j,1,e) &
 
  966                        + dx(i,3) * u(3,j,1,e) &
 
  967                        + dx(i,4) * u(4,j,1,e) &
 
  968                        + dx(i,5) * u(5,j,1,e) &
 
  969                        + dx(i,6) * u(6,j,1,e) &
 
  970                        + dx(i,7) * u(7,j,1,e) &
 
  971                        + dx(i,8) * u(8,j,1,e) &
 
  972                        + dx(i,9) * u(9,j,1,e) &
 
  973                        + dx(i,10) * u(10,j,1,e) &
 
  974                        + dx(i,11) * u(11,j,1,e)
 
  981                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
  982                          + dy(j,2) * u(i,2,k,e) &
 
  983                          + dy(j,3) * u(i,3,k,e) &
 
  984                          + dy(j,4) * u(i,4,k,e) &
 
  985                          + dy(j,5) * u(i,5,k,e) &
 
  986                          + dy(j,6) * u(i,6,k,e) &
 
  987                          + dy(j,7) * u(i,7,k,e) &
 
  988                          + dy(j,8) * u(i,8,k,e) &
 
  989                          + dy(j,9) * u(i,9,k,e) &
 
  990                          + dy(j,10) * u(i,10,k,e) &
 
  991                          + dy(j,11) * u(i,11,k,e)
 
  998             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
  999                       + dz(k,2) * u(i,1,2,e) &
 
 1000                       + dz(k,3) * u(i,1,3,e) &
 
 1001                       + dz(k,4) * u(i,1,4,e) &
 
 1002                       + dz(k,5) * u(i,1,5,e) &
 
 1003                       + dz(k,6) * u(i,1,6,e) &
 
 1004                       + dz(k,7) * u(i,1,7,e) &
 
 1005                       + dz(k,8) * u(i,1,8,e) &
 
 1006                       + dz(k,9) * u(i,1,9,e) &
 
 1007                       + dz(k,10) * u(i,1,10,e) &
 
 1008                       + dz(k,11) * u(i,1,11,e)
 
 1012       do i = 1, lx * lx * lx
 
 1013          ux(i,1,1,e) = w3(i,1,1) &
 
 1014                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
 1015                        + dsdx(i,1,1,e) * us(i,1,1) &
 
 1016                        + dtdx(i,1,1,e) * ut(i,1,1) )
 
 1017          uy(i,1,1,e) = w3(i,1,1) &
 
 1018                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
 1019                        + drdy(i,1,1,e) * ur(i,1,1) &
 
 1020                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
 1021          uz(i,1,1,e) = w3(i,1,1) &
 
 1022                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
 1023                        + drdz(i,1,1,e) * ur(i,1,1) &
 
 1024                        + dsdz(i,1,1,e) * us(i,1,1) )
 
 1027  end subroutine cpu_opgrad_lx11
 
 1029  subroutine cpu_opgrad_lx10(ux, uy, uz, u, dx, dy, dz, &
 
 1030       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
 1031    integer, 
parameter :: lx = 10
 
 1032    integer, 
intent(in) :: n
 
 1033    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
 1034    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
 1035    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 1036    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
 1037    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
 1038    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
 1039    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 1040    real(kind=rp) :: ur(lx, lx, lx)
 
 1041    real(kind=rp) :: us(lx, lx, lx)
 
 1042    real(kind=rp) :: ut(lx, lx, lx)
 
 1043    integer :: e, i, j, k
 
 1048             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
 1049                        + dx(i,2) * u(2,j,1,e) &
 
 1050                        + dx(i,3) * u(3,j,1,e) &
 
 1051                        + dx(i,4) * u(4,j,1,e) &
 
 1052                        + dx(i,5) * u(5,j,1,e) &
 
 1053                        + dx(i,6) * u(6,j,1,e) &
 
 1054                        + dx(i,7) * u(7,j,1,e) &
 
 1055                        + dx(i,8) * u(8,j,1,e) &
 
 1056                        + dx(i,9) * u(9,j,1,e) &
 
 1057                        + dx(i,10) * u(10,j,1,e)
 
 1064                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
 1065                          + dy(j,2) * u(i,2,k,e) &
 
 1066                          + dy(j,3) * u(i,3,k,e) &
 
 1067                          + dy(j,4) * u(i,4,k,e) &
 
 1068                          + dy(j,5) * u(i,5,k,e) &
 
 1069                          + dy(j,6) * u(i,6,k,e) &
 
 1070                          + dy(j,7) * u(i,7,k,e) &
 
 1071                          + dy(j,8) * u(i,8,k,e) &
 
 1072                          + dy(j,9) * u(i,9,k,e) &
 
 1073                          + dy(j,10) * u(i,10,k,e)
 
 1080             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
 1081                       + dz(k,2) * u(i,1,2,e) &
 
 1082                       + dz(k,3) * u(i,1,3,e) &
 
 1083                       + dz(k,4) * u(i,1,4,e) &
 
 1084                       + dz(k,5) * u(i,1,5,e) &
 
 1085                       + dz(k,6) * u(i,1,6,e) &
 
 1086                       + dz(k,7) * u(i,1,7,e) &
 
 1087                       + dz(k,8) * u(i,1,8,e) &
 
 1088                       + dz(k,9) * u(i,1,9,e) &
 
 1089                       + dz(k,10) * u(i,1,10,e)
 
 1093       do i = 1, lx * lx * lx
 
 1094          ux(i,1,1,e) = w3(i,1,1) &
 
 1095                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
 1096                         + dsdx(i,1,1,e) * us(i,1,1) &
 
 1097                         + dtdx(i,1,1,e) * ut(i,1,1) )
 
 1098          uy(i,1,1,e) = w3(i,1,1) &
 
 1099                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
 1100                        + drdy(i,1,1,e) * ur(i,1,1) &
 
 1101                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
 1102          uz(i,1,1,e) = w3(i,1,1) &
 
 1103                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
 1104                        + drdz(i,1,1,e) * ur(i,1,1) &
 
 1105                        + dsdz(i,1,1,e) * us(i,1,1) )
 
 1108  end subroutine cpu_opgrad_lx10
 
 1110  subroutine cpu_opgrad_lx9(ux, uy, uz, u, dx, dy, dz, &
 
 1111       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
 1112    integer, 
parameter :: lx = 9
 
 1113    integer, 
intent(in) :: n
 
 1114    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
 1115    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
 1116    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 1117    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
 1118    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
 1119    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
 1120    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 1121    real(kind=rp) :: ur(lx, lx, lx)
 
 1122    real(kind=rp) :: us(lx, lx, lx)
 
 1123    real(kind=rp) :: ut(lx, lx, lx)
 
 1124    integer :: e, i, j, k
 
 1129             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
 1130                        + dx(i,2) * u(2,j,1,e) &
 
 1131                        + dx(i,3) * u(3,j,1,e) &
 
 1132                        + dx(i,4) * u(4,j,1,e) &
 
 1133                        + dx(i,5) * u(5,j,1,e) &
 
 1134                        + dx(i,6) * u(6,j,1,e) &
 
 1135                        + dx(i,7) * u(7,j,1,e) &
 
 1136                        + dx(i,8) * u(8,j,1,e) &
 
 1137                        + dx(i,9) * u(9,j,1,e)
 
 1144                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
 1145                          + dy(j,2) * u(i,2,k,e) &
 
 1146                          + dy(j,3) * u(i,3,k,e) &
 
 1147                          + dy(j,4) * u(i,4,k,e) &
 
 1148                          + dy(j,5) * u(i,5,k,e) &
 
 1149                          + dy(j,6) * u(i,6,k,e) &
 
 1150                          + dy(j,7) * u(i,7,k,e) &
 
 1151                          + dy(j,8) * u(i,8,k,e) &
 
 1152                          + dy(j,9) * u(i,9,k,e)
 
 1159             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
 1160                       + dz(k,2) * u(i,1,2,e) &
 
 1161                       + dz(k,3) * u(i,1,3,e) &
 
 1162                       + dz(k,4) * u(i,1,4,e) &
 
 1163                       + dz(k,5) * u(i,1,5,e) &
 
 1164                       + dz(k,6) * u(i,1,6,e) &
 
 1165                       + dz(k,7) * u(i,1,7,e) &
 
 1166                       + dz(k,8) * u(i,1,8,e) &
 
 1167                       + dz(k,9) * u(i,1,9,e)
 
 1171       do i = 1, lx * lx * lx
 
 1172          ux(i,1,1,e) = w3(i,1,1) &
 
 1173                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
 1174                        + dsdx(i,1,1,e) * us(i,1,1) &
 
 1175                        + dtdx(i,1,1,e) * ut(i,1,1) )
 
 1176          uy(i,1,1,e) = w3(i,1,1) &
 
 1177                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
 1178                        + drdy(i,1,1,e) * ur(i,1,1) &
 
 1179                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
 1180          uz(i,1,1,e) = w3(i,1,1) &
 
 1181                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
 1182                        + drdz(i,1,1,e) * ur(i,1,1) &
 
 1183                        + dsdz(i,1,1,e) * us(i,1,1) )
 
 1186  end subroutine cpu_opgrad_lx9
 
 1188  subroutine cpu_opgrad_lx8(ux, uy, uz, u, dx, dy, dz, &
 
 1189       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
 1190    integer, 
parameter :: lx = 8
 
 1191    integer, 
intent(in) :: n
 
 1192    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
 1193    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
 1194    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 1195    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
 1196    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
 1197    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
 1198    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 1199    real(kind=rp) :: ur(lx, lx, lx)
 
 1200    real(kind=rp) :: us(lx, lx, lx)
 
 1201    real(kind=rp) :: ut(lx, lx, lx)
 
 1202    integer :: e, i, j, k
 
 1207             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
 1208                        + dx(i,2) * u(2,j,1,e) &
 
 1209                        + dx(i,3) * u(3,j,1,e) &
 
 1210                        + dx(i,4) * u(4,j,1,e) &
 
 1211                        + dx(i,5) * u(5,j,1,e) &
 
 1212                        + dx(i,6) * u(6,j,1,e) &
 
 1213                        + dx(i,7) * u(7,j,1,e) &
 
 1214                        + dx(i,8) * u(8,j,1,e)
 
 1221                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
 1222                          + dy(j,2) * u(i,2,k,e) &
 
 1223                          + dy(j,3) * u(i,3,k,e) &
 
 1224                          + dy(j,4) * u(i,4,k,e) &
 
 1225                          + dy(j,5) * u(i,5,k,e) &
 
 1226                          + dy(j,6) * u(i,6,k,e) &
 
 1227                          + dy(j,7) * u(i,7,k,e) &
 
 1228                          + dy(j,8) * u(i,8,k,e)
 
 1235             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
 1236                       + dz(k,2) * u(i,1,2,e) &
 
 1237                       + dz(k,3) * u(i,1,3,e) &
 
 1238                       + dz(k,4) * u(i,1,4,e) &
 
 1239                       + dz(k,5) * u(i,1,5,e) &
 
 1240                       + dz(k,6) * u(i,1,6,e) &
 
 1241                       + dz(k,7) * u(i,1,7,e) &
 
 1242                       + dz(k,8) * u(i,1,8,e)
 
 1246       do i = 1, lx * lx * lx
 
 1247          ux(i,1,1,e) = w3(i,1,1) &
 
 1248                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
 1249                        + dsdx(i,1,1,e) * us(i,1,1) &
 
 1250                        + dtdx(i,1,1,e) * ut(i,1,1) )
 
 1251          uy(i,1,1,e) = w3(i,1,1) &
 
 1252                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
 1253                        + drdy(i,1,1,e) * ur(i,1,1) &
 
 1254                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
 1255          uz(i,1,1,e) = w3(i,1,1) &
 
 1256                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
 1257                        + drdz(i,1,1,e) * ur(i,1,1) &
 
 1258                        + dsdz(i,1,1,e) * us(i,1,1) )
 
 1261  end subroutine cpu_opgrad_lx8
 
 1263  subroutine cpu_opgrad_lx7(ux, uy, uz, u, dx, dy, dz, &
 
 1264       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
 1265    integer, 
parameter :: lx = 7
 
 1266    integer, 
intent(in) :: n
 
 1267    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
 1268    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
 1269    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 1270    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
 1271    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
 1272    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
 1273    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 1274    real(kind=rp) :: ur(lx, lx, lx)
 
 1275    real(kind=rp) :: us(lx, lx, lx)
 
 1276    real(kind=rp) :: ut(lx, lx, lx)
 
 1277    integer :: e, i, j, k
 
 1282             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
 1283                        + dx(i,2) * u(2,j,1,e) &
 
 1284                        + dx(i,3) * u(3,j,1,e) &
 
 1285                        + dx(i,4) * u(4,j,1,e) &
 
 1286                        + dx(i,5) * u(5,j,1,e) &
 
 1287                        + dx(i,6) * u(6,j,1,e) &
 
 1288                        + dx(i,7) * u(7,j,1,e)
 
 1295                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
 1296                          + dy(j,2) * u(i,2,k,e) &
 
 1297                          + dy(j,3) * u(i,3,k,e) &
 
 1298                          + dy(j,4) * u(i,4,k,e) &
 
 1299                          + dy(j,5) * u(i,5,k,e) &
 
 1300                          + dy(j,6) * u(i,6,k,e) &
 
 1301                          + dy(j,7) * u(i,7,k,e)
 
 1308             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
 1309                       + dz(k,2) * u(i,1,2,e) &
 
 1310                       + dz(k,3) * u(i,1,3,e) &
 
 1311                       + dz(k,4) * u(i,1,4,e) &
 
 1312                       + dz(k,5) * u(i,1,5,e) &
 
 1313                       + dz(k,6) * u(i,1,6,e) &
 
 1314                       + dz(k,7) * u(i,1,7,e)
 
 1318       do i = 1, lx * lx * lx
 
 1319          ux(i,1,1,e) = w3(i,1,1) &
 
 1320                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
 1321                        + dsdx(i,1,1,e) * us(i,1,1) &
 
 1322                        + dtdx(i,1,1,e) * ut(i,1,1) )
 
 1323          uy(i,1,1,e) = w3(i,1,1) &
 
 1324                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
 1325                        + drdy(i,1,1,e) * ur(i,1,1) &
 
 1326                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
 1327          uz(i,1,1,e) = w3(i,1,1) &
 
 1328                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
 1329                        + drdz(i,1,1,e) * ur(i,1,1) &
 
 1330                        + dsdz(i,1,1,e) * us(i,1,1) )
 
 1333  end subroutine cpu_opgrad_lx7
 
 1335  subroutine cpu_opgrad_lx6(ux, uy, uz, u, dx, dy, dz, &
 
 1336       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
 1337    integer, 
parameter :: lx = 6
 
 1338    integer, 
intent(in) :: n
 
 1339    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
 1340    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
 1341    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 1342    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
 1343    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
 1344    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
 1345    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 1346    real(kind=rp) :: ur(lx, lx, lx)
 
 1347    real(kind=rp) :: us(lx, lx, lx)
 
 1348    real(kind=rp) :: ut(lx, lx, lx)
 
 1349    integer :: e, i, j, k
 
 1354             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
 1355                        + dx(i,2) * u(2,j,1,e) &
 
 1356                        + dx(i,3) * u(3,j,1,e) &
 
 1357                        + dx(i,4) * u(4,j,1,e) &
 
 1358                        + dx(i,5) * u(5,j,1,e) &
 
 1359                        + dx(i,6) * u(6,j,1,e)
 
 1366                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
 1367                          + dy(j,2) * u(i,2,k,e) &
 
 1368                          + dy(j,3) * u(i,3,k,e) &
 
 1369                          + dy(j,4) * u(i,4,k,e) &
 
 1370                          + dy(j,5) * u(i,5,k,e) &
 
 1371                          + dy(j,6) * u(i,6,k,e)
 
 1378             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
 1379                       + dz(k,2) * u(i,1,2,e) &
 
 1380                       + dz(k,3) * u(i,1,3,e) &
 
 1381                       + dz(k,4) * u(i,1,4,e) &
 
 1382                       + dz(k,5) * u(i,1,5,e) &
 
 1383                       + dz(k,6) * u(i,1,6,e)
 
 1387       do i = 1, lx * lx * lx
 
 1388          ux(i,1,1,e) = w3(i,1,1) &
 
 1389                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
 1390                        + dsdx(i,1,1,e) * us(i,1,1) &
 
 1391                        + dtdx(i,1,1,e) * ut(i,1,1) )
 
 1392          uy(i,1,1,e) = w3(i,1,1) &
 
 1393                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
 1394                        + drdy(i,1,1,e) * ur(i,1,1) &
 
 1395                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
 1396          uz(i,1,1,e) = w3(i,1,1) &
 
 1397                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
 1398                        + drdz(i,1,1,e) * ur(i,1,1) &
 
 1399                        + dsdz(i,1,1,e) * us(i,1,1) )
 
 1402  end subroutine cpu_opgrad_lx6
 
 1404  subroutine cpu_opgrad_lx5(ux, uy, uz, u, dx, dy, dz, &
 
 1405       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
 1406    integer, 
parameter :: lx = 5
 
 1407    integer, 
intent(in) :: n
 
 1408    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
 1409    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
 1410    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 1411    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
 1412    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
 1413    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
 1414    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 1415    real(kind=rp) :: ur(lx, lx, lx)
 
 1416    real(kind=rp) :: us(lx, lx, lx)
 
 1417    real(kind=rp) :: ut(lx, lx, lx)
 
 1418    integer :: e, i, j, k
 
 1423             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
 1424                        + dx(i,2) * u(2,j,1,e) &
 
 1425                        + dx(i,3) * u(3,j,1,e) &
 
 1426                        + dx(i,4) * u(4,j,1,e) &
 
 1427                        + dx(i,5) * u(5,j,1,e)
 
 1434                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
 1435                          + dy(j,2) * u(i,2,k,e) &
 
 1436                          + dy(j,3) * u(i,3,k,e) &
 
 1437                          + dy(j,4) * u(i,4,k,e) &
 
 1438                          + dy(j,5) * u(i,5,k,e)
 
 1445             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
 1446                       + dz(k,2) * u(i,1,2,e) &
 
 1447                       + dz(k,3) * u(i,1,3,e) &
 
 1448                       + dz(k,4) * u(i,1,4,e) &
 
 1449                       + dz(k,5) * u(i,1,5,e)
 
 1453       do i = 1, lx * lx * lx
 
 1454          ux(i,1,1,e) = w3(i,1,1) &
 
 1455                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
 1456                        + dsdx(i,1,1,e) * us(i,1,1) &
 
 1457                        + dtdx(i,1,1,e) * ut(i,1,1) )
 
 1458          uy(i,1,1,e) = w3(i,1,1) &
 
 1459                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
 1460                        + drdy(i,1,1,e) * ur(i,1,1) &
 
 1461                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
 1462          uz(i,1,1,e) = w3(i,1,1) &
 
 1463                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
 1464                        + drdz(i,1,1,e) * ur(i,1,1) &
 
 1465                        + dsdz(i,1,1,e) * us(i,1,1) )
 
 1468  end subroutine cpu_opgrad_lx5
 
 1470  subroutine cpu_opgrad_lx4(ux, uy, uz, u, dx, dy, dz, &
 
 1471       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
 1472    integer, 
parameter :: lx = 4
 
 1473    integer, 
intent(in) :: n
 
 1474    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
 1475    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
 1476    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 1477    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
 1478    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
 1479    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
 1480    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 1481    real(kind=rp) :: ur(lx, lx, lx)
 
 1482    real(kind=rp) :: us(lx, lx, lx)
 
 1483    real(kind=rp) :: ut(lx, lx, lx)
 
 1484    integer :: e, i, j, k
 
 1489             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
 1490                        + dx(i,2) * u(2,j,1,e) &
 
 1491                        + dx(i,3) * u(3,j,1,e) &
 
 1492                        + dx(i,4) * u(4,j,1,e)
 
 1499                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
 1500                          + dy(j,2) * u(i,2,k,e) &
 
 1501                          + dy(j,3) * u(i,3,k,e) &
 
 1502                          + dy(j,4) * u(i,4,k,e)
 
 1509             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
 1510                       + dz(k,2) * u(i,1,2,e) &
 
 1511                       + dz(k,3) * u(i,1,3,e) &
 
 1512                       + dz(k,4) * u(i,1,4,e)
 
 1516       do i = 1, lx * lx * lx
 
 1517          ux(i,1,1,e) = w3(i,1,1) &
 
 1518                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
 1519                        + dsdx(i,1,1,e) * us(i,1,1) &
 
 1520                        + dtdx(i,1,1,e) * ut(i,1,1) )
 
 1521          uy(i,1,1,e) = w3(i,1,1) &
 
 1522                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
 1523                        + drdy(i,1,1,e) * ur(i,1,1) &
 
 1524                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
 1525          uz(i,1,1,e) = w3(i,1,1) &
 
 1526                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
 1527                        + drdz(i,1,1,e) * ur(i,1,1) &
 
 1528                        + dsdz(i,1,1,e) * us(i,1,1) )
 
 1531  end subroutine cpu_opgrad_lx4
 
 1533  subroutine cpu_opgrad_lx3(ux, uy, uz, u, dx, dy, dz, &
 
 1534       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
 1535    integer, 
parameter :: lx = 3
 
 1536    integer, 
intent(in) :: n
 
 1537    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
 1538    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
 1539    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 1540    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
 1541    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
 1542    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
 1543    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 1544    real(kind=rp) :: ur(lx, lx, lx)
 
 1545    real(kind=rp) :: us(lx, lx, lx)
 
 1546    real(kind=rp) :: ut(lx, lx, lx)
 
 1547    integer :: e, i, j, k
 
 1552             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
 1553                        + dx(i,2) * u(2,j,1,e) &
 
 1554                        + dx(i,3) * u(3,j,1,e)
 
 1561                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
 1562                          + dy(j,2) * u(i,2,k,e) &
 
 1563                          + dy(j,3) * u(i,3,k,e)
 
 1570             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
 1571                       + dz(k,2) * u(i,1,2,e) &
 
 1572                       + dz(k,3) * u(i,1,3,e)
 
 1576       do i = 1, lx * lx * lx
 
 1577          ux(i,1,1,e) = w3(i,1,1) &
 
 1578                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
 1579                        + dsdx(i,1,1,e) * us(i,1,1) &
 
 1580                        + dtdx(i,1,1,e) * ut(i,1,1) )
 
 1581          uy(i,1,1,e) = w3(i,1,1) &
 
 1582                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
 1583                        + drdy(i,1,1,e) * ur(i,1,1) &
 
 1584                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
 1585          uz(i,1,1,e) = w3(i,1,1) &
 
 1586                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
 1587                        + drdz(i,1,1,e) * ur(i,1,1) &
 
 1588                        + dsdz(i,1,1,e) * us(i,1,1) )
 
 1591  end subroutine cpu_opgrad_lx3
 
 1593  subroutine cpu_opgrad_lx2(ux, uy, uz, u, dx, dy, dz, &
 
 1594       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
 
 1595    integer, 
parameter :: lx = 2
 
 1596    integer, 
intent(in) :: n
 
 1597    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(inout) :: ux, uy, uz
 
 1598    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: u
 
 1599    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 1600    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdx, dsdx, dtdx
 
 1601    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdy, dsdy, dtdy
 
 1602    real(kind=rp), 
dimension(lx, lx, lx, n), 
intent(in) :: drdz, dsdz, dtdz
 
 1603    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 1604    real(kind=rp) :: ur(lx, lx, lx)
 
 1605    real(kind=rp) :: us(lx, lx, lx)
 
 1606    real(kind=rp) :: ut(lx, lx, lx)
 
 1607    integer :: e, i, j, k
 
 1612             ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
 
 1613                        + dx(i,2) * u(2,j,1,e)
 
 1620                us(i,j,k) = dy(j,1) * u(i,1,k,e) &
 
 1621                          + dy(j,2) * u(i,2,k,e)
 
 1628             ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
 
 1629                       + dz(k,2) * u(i,1,2,e)
 
 1633       do i = 1, lx * lx * lx
 
 1634          ux(i,1,1,e) = w3(i,1,1) &
 
 1635                      * ( drdx(i,1,1,e) * ur(i,1,1) &
 
 1636                        + dsdx(i,1,1,e) * us(i,1,1) &
 
 1637                        + dtdx(i,1,1,e) * ut(i,1,1) )
 
 1638          uy(i,1,1,e) = w3(i,1,1) &
 
 1639                      * ( dsdy(i,1,1,e) * us(i,1,1) &
 
 1640                        + drdy(i,1,1,e) * ur(i,1,1) &
 
 1641                        + dtdy(i,1,1,e) * ut(i,1,1) )
 
 1642          uz(i,1,1,e) = w3(i,1,1) &
 
 1643                      * ( dtdz(i,1,1,e) * ut(i,1,1) &
 
 1644                        + drdz(i,1,1,e) * ur(i,1,1) &
 
 1645                        + dsdz(i,1,1,e) * us(i,1,1) )
 
 1648  end subroutine cpu_opgrad_lx2
 
 1650  subroutine opr_cpu_opgrad_single(ux, uy, uz, u, coef, e)
 
 1651    integer, 
parameter :: e_len = 1
 
 1652    type(coef_t), 
intent(in) :: coef
 
 1653    integer, 
intent(in) :: e
 
 1654    real(kind=rp), 
dimension(coef%Xh%lxyz, e_len), 
intent(inout) :: ux
 
 1655    real(kind=rp), 
dimension(coef%Xh%lxyz, e_len), 
intent(inout) :: uy
 
 1656    real(kind=rp), 
dimension(coef%Xh%lxyz, e_len), 
intent(inout) :: uz
 
 1657    real(kind=rp), 
dimension(coef%Xh%lxyz, e_len), 
intent(in) :: u
 
 1659    associate(xh => coef%Xh, msh => coef%msh, &
 
 1660         drdx => coef%drdx, drdy => coef%drdy, drdz => coef%drdz, &
 
 1661         dsdx => coef%dsdx, dsdy => coef%dsdy, dsdz => coef%dsdz, &
 
 1662         dtdx => coef%dtdx, dtdy => coef%dtdy, dtdz => coef%dtdz)
 
 1666         call cpu_opgrad_lx18_single(ux, uy, uz, u, &
 
 1667              xh%dx, xh%dy, xh%dz, &
 
 1668              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1669              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1670              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1673         call cpu_opgrad_lx17_single(ux, uy, uz, u, &
 
 1674              xh%dx, xh%dy, xh%dz, &
 
 1675              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1676              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1677              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1680         call cpu_opgrad_lx16_single(ux, uy, uz, u, &
 
 1681              xh%dx, xh%dy, xh%dz, &
 
 1682              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1683              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1684              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1687         call cpu_opgrad_lx15_single(ux, uy, uz, u, &
 
 1688              xh%dx, xh%dy, xh%dz, &
 
 1689              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1690              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1691              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1694         call cpu_opgrad_lx14_single(ux, uy, uz, u, &
 
 1695              xh%dx, xh%dy, xh%dz, &
 
 1696              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1697              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1698              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1701         call cpu_opgrad_lx13_single(ux, uy, uz, u, &
 
 1702              xh%dx, xh%dy, xh%dz, &
 
 1703              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1704              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1705              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1708         call cpu_opgrad_lx12_single(ux, uy, uz, u, &
 
 1709              xh%dx, xh%dy, xh%dz, &
 
 1710              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1711              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1712              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1715         call cpu_opgrad_lx11_single(ux, uy, uz, u, &
 
 1716              xh%dx, xh%dy, xh%dz, &
 
 1717              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1718              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1719              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1722         call cpu_opgrad_lx10_single(ux, uy, uz, u, &
 
 1723              xh%dx, xh%dy, xh%dz, &
 
 1724              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1725              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1726              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1730         call cpu_opgrad_lx9_single(ux, uy, uz, u, &
 
 1731              xh%dx, xh%dy, xh%dz, &
 
 1732              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1733              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1734              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1737         call cpu_opgrad_lx8_single(ux, uy, uz, u, &
 
 1738              xh%dx, xh%dy, xh%dz, &
 
 1739              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1740              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1741              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1744         call cpu_opgrad_lx7_single(ux, uy, uz, u, &
 
 1745              xh%dx, xh%dy, xh%dz, &
 
 1746              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1747              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1748              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1751         call cpu_opgrad_lx6_single(ux, uy, uz, u, &
 
 1752              xh%dx, xh%dy, xh%dz, &
 
 1753              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1754              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1755              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1758         call cpu_opgrad_lx5_single(ux, uy, uz, u, &
 
 1759              xh%dx, xh%dy, xh%dz, &
 
 1760              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1761              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1762              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1765         call cpu_opgrad_lx4_single(ux, uy, uz, u, &
 
 1766              xh%dx, xh%dy, xh%dz, &
 
 1767              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1768              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1769              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1772         call cpu_opgrad_lx3_single(ux, uy, uz, u, &
 
 1773              xh%dx, xh%dy, xh%dz, &
 
 1774              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1775              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1776              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1779         call cpu_opgrad_lx2_single(ux, uy, uz, u, &
 
 1780              xh%dx, xh%dy, xh%dz, &
 
 1781              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1782              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1783              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1786         call cpu_opgrad_lx_single(ux, uy, uz, u, &
 
 1787              xh%dx, xh%dy, xh%dz, &
 
 1788              drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
 
 1789              drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
 
 1790              drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
 
 1795  end subroutine opr_cpu_opgrad_single
 
 1797  subroutine cpu_opgrad_lx_single(ux, uy, uz, u, dx, dy, dz, &
 
 1798       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, lx)
 
 1799    integer, 
intent(in) :: lx
 
 1800    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 1801    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 1802    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 1803    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 1804    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 1805    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 1806    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 1807    real(kind=rp) :: ur(lx, lx, lx)
 
 1808    real(kind=rp) :: us(lx, lx, lx)
 
 1809    real(kind=rp) :: ut(lx, lx, lx)
 
 1810    real(kind=rp) :: tmp
 
 1811    integer :: i, j, k, l
 
 1817             tmp = tmp + dx(i,k) * u(k,j,1)
 
 1828                tmp = tmp + dy(j,l) * u(i,l,k)
 
 1839             tmp = tmp + dz(k,l) * u(i,1,l)
 
 1845    do i = 1, lx * lx * lx
 
 1846       ux(i,1,1) = w3(i,1,1) &
 
 1847                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 1848                     + dsdx(i,1,1) * us(i,1,1) &
 
 1849                     + dtdx(i,1,1) * ut(i,1,1) )
 
 1850       uy(i,1,1) = w3(i,1,1) &
 
 1851                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 1852                     + drdy(i,1,1) * ur(i,1,1) &
 
 1853                     + dtdy(i,1,1) * ut(i,1,1) )
 
 1854       uz(i,1,1) = w3(i,1,1) &
 
 1855                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 1856                     + drdz(i,1,1) * ur(i,1,1) &
 
 1857                     + dsdz(i,1,1) * us(i,1,1) )
 
 1860  end subroutine cpu_opgrad_lx_single
 
 1862  subroutine cpu_opgrad_lx18_single(ux, uy, uz, u, dx, dy, dz, &
 
 1863       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 1864    integer, 
parameter :: lx = 18
 
 1865      real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 1866    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 1867    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 1868    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 1869    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 1870    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 1871    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 1872    real(kind=rp) :: ur(lx, lx, lx)
 
 1873    real(kind=rp) :: us(lx, lx, lx)
 
 1874    real(kind=rp) :: ut(lx, lx, lx)
 
 1879          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 1880                    + dx(i,2) * u(2,j,1) &
 
 1881                    + dx(i,3) * u(3,j,1) &
 
 1882                    + dx(i,4) * u(4,j,1) &
 
 1883                    + dx(i,5) * u(5,j,1) &
 
 1884                    + dx(i,6) * u(6,j,1) &
 
 1885                    + dx(i,7) * u(7,j,1) &
 
 1886                    + dx(i,8) * u(8,j,1) &
 
 1887                    + dx(i,9) * u(9,j,1) &
 
 1888                    + dx(i,10) * u(10,j,1) &
 
 1889                    + dx(i,11) * u(11,j,1) &
 
 1890                    + dx(i,12) * u(12,j,1) &
 
 1891                    + dx(i,13) * u(13,j,1) &
 
 1892                    + dx(i,14) * u(14,j,1) &
 
 1893                    + dx(i,15) * u(15,j,1) &
 
 1894                    + dx(i,16) * u(16,j,1) &
 
 1895                    + dx(i,17) * u(17,j,1) &
 
 1896                    + dx(i,18) * u(18,j,1)
 
 1903             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 1904                       + dy(j,2) * u(i,2,k) &
 
 1905                       + dy(j,3) * u(i,3,k) &
 
 1906                       + dy(j,4) * u(i,4,k) &
 
 1907                       + dy(j,5) * u(i,5,k) &
 
 1908                       + dy(j,6) * u(i,6,k) &
 
 1909                       + dy(j,7) * u(i,7,k) &
 
 1910                       + dy(j,8) * u(i,8,k) &
 
 1911                       + dy(j,9) * u(i,9,k) &
 
 1912                       + dy(j,10) * u(i,10,k) &
 
 1913                       + dy(j,11) * u(i,11,k) &
 
 1914                       + dy(j,12) * u(i,12,k) &
 
 1915                       + dy(j,13) * u(i,13,k) &
 
 1916                       + dy(j,14) * u(i,14,k) &
 
 1917                       + dy(j,15) * u(i,15,k) &
 
 1918                       + dy(j,16) * u(i,16,k) &
 
 1919                       + dy(j,17) * u(i,17,k) &
 
 1920                       + dy(j,18) * u(i,18,k)
 
 1927          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 1928                    + dz(k,2) * u(i,1,2) &
 
 1929                    + dz(k,3) * u(i,1,3) &
 
 1930                    + dz(k,4) * u(i,1,4) &
 
 1931                    + dz(k,5) * u(i,1,5) &
 
 1932                    + dz(k,6) * u(i,1,6) &
 
 1933                    + dz(k,7) * u(i,1,7) &
 
 1934                    + dz(k,8) * u(i,1,8) &
 
 1935                    + dz(k,9) * u(i,1,9) &
 
 1936                    + dz(k,10) * u(i,1,10) &
 
 1937                    + dz(k,11) * u(i,1,11) &
 
 1938                    + dz(k,12) * u(i,1,12) &
 
 1939                    + dz(k,13) * u(i,1,13) &
 
 1940                    + dz(k,14) * u(i,1,14) &
 
 1941                    + dz(k,15) * u(i,1,15) &
 
 1942                    + dz(k,16) * u(i,1,16) &
 
 1943                    + dz(k,17) * u(i,1,17) &
 
 1944                    + dz(k,18) * u(i,1,18)
 
 1948    do i = 1, lx * lx * lx
 
 1949       ux(i,1,1) = w3(i,1,1) &
 
 1950                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 1951                     + dsdx(i,1,1) * us(i,1,1) &
 
 1952                     + dtdx(i,1,1) * ut(i,1,1) )
 
 1953       uy(i,1,1) = w3(i,1,1) &
 
 1954                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 1955                     + drdy(i,1,1) * ur(i,1,1) &
 
 1956                     + dtdy(i,1,1) * ut(i,1,1) )
 
 1957       uz(i,1,1) = w3(i,1,1) &
 
 1958                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 1959                     + drdz(i,1,1) * ur(i,1,1) &
 
 1960                     + dsdz(i,1,1) * us(i,1,1) )
 
 1963  end subroutine cpu_opgrad_lx18_single
 
 1965  subroutine cpu_opgrad_lx17_single(ux, uy, uz, u, dx, dy, dz, &
 
 1966       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 1967    integer, 
parameter :: lx = 17
 
 1968    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 1969    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 1970    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 1971    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 1972    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 1973    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 1974    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 1975    real(kind=rp) :: ur(lx, lx, lx)
 
 1976    real(kind=rp) :: us(lx, lx, lx)
 
 1977    real(kind=rp) :: ut(lx, lx, lx)
 
 1982          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 1983                    + dx(i,2) * u(2,j,1) &
 
 1984                    + dx(i,3) * u(3,j,1) &
 
 1985                    + dx(i,4) * u(4,j,1) &
 
 1986                    + dx(i,5) * u(5,j,1) &
 
 1987                    + dx(i,6) * u(6,j,1) &
 
 1988                    + dx(i,7) * u(7,j,1) &
 
 1989                    + dx(i,8) * u(8,j,1) &
 
 1990                    + dx(i,9) * u(9,j,1) &
 
 1991                    + dx(i,10) * u(10,j,1) &
 
 1992                    + dx(i,11) * u(11,j,1) &
 
 1993                    + dx(i,12) * u(12,j,1) &
 
 1994                    + dx(i,13) * u(13,j,1) &
 
 1995                    + dx(i,14) * u(14,j,1) &
 
 1996                    + dx(i,15) * u(15,j,1) &
 
 1997                    + dx(i,16) * u(16,j,1) &
 
 1998                    + dx(i,17) * u(17,j,1)
 
 2005             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2006                       + dy(j,2) * u(i,2,k) &
 
 2007                       + dy(j,3) * u(i,3,k) &
 
 2008                       + dy(j,4) * u(i,4,k) &
 
 2009                       + dy(j,5) * u(i,5,k) &
 
 2010                       + dy(j,6) * u(i,6,k) &
 
 2011                       + dy(j,7) * u(i,7,k) &
 
 2012                       + dy(j,8) * u(i,8,k) &
 
 2013                       + dy(j,9) * u(i,9,k) &
 
 2014                       + dy(j,10) * u(i,10,k) &
 
 2015                       + dy(j,11) * u(i,11,k) &
 
 2016                       + dy(j,12) * u(i,12,k) &
 
 2017                       + dy(j,13) * u(i,13,k) &
 
 2018                       + dy(j,14) * u(i,14,k) &
 
 2019                       + dy(j,15) * u(i,15,k) &
 
 2020                       + dy(j,16) * u(i,16,k) &
 
 2021                       + dy(j,17) * u(i,17,k)
 
 2028          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2029                    + dz(k,2) * u(i,1,2) &
 
 2030                    + dz(k,3) * u(i,1,3) &
 
 2031                    + dz(k,4) * u(i,1,4) &
 
 2032                    + dz(k,5) * u(i,1,5) &
 
 2033                    + dz(k,6) * u(i,1,6) &
 
 2034                    + dz(k,7) * u(i,1,7) &
 
 2035                    + dz(k,8) * u(i,1,8) &
 
 2036                    + dz(k,9) * u(i,1,9) &
 
 2037                    + dz(k,10) * u(i,1,10) &
 
 2038                    + dz(k,11) * u(i,1,11) &
 
 2039                    + dz(k,12) * u(i,1,12) &
 
 2040                    + dz(k,13) * u(i,1,13) &
 
 2041                    + dz(k,14) * u(i,1,14) &
 
 2042                    + dz(k,15) * u(i,1,15) &
 
 2043                    + dz(k,16) * u(i,1,16) &
 
 2044                    + dz(k,17) * u(i,1,17)
 
 2048    do i = 1, lx * lx * lx
 
 2049       ux(i,1,1) = w3(i,1,1) &
 
 2050                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 2051                     + dsdx(i,1,1) * us(i,1,1) &
 
 2052                     + dtdx(i,1,1) * ut(i,1,1) )
 
 2053       uy(i,1,1) = w3(i,1,1) &
 
 2054                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 2055                     + drdy(i,1,1) * ur(i,1,1) &
 
 2056                     + dtdy(i,1,1) * ut(i,1,1) )
 
 2057       uz(i,1,1) = w3(i,1,1) &
 
 2058                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 2059                     + drdz(i,1,1) * ur(i,1,1) &
 
 2060                     + dsdz(i,1,1) * us(i,1,1) )
 
 2062  end subroutine cpu_opgrad_lx17_single
 
 2064  subroutine cpu_opgrad_lx16_single(ux, uy, uz, u, dx, dy, dz, &
 
 2065       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 2066    integer, 
parameter :: lx = 16
 
 2067    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 2068    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 2069    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 2070    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 2071    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 2072    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 2073    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 2074    real(kind=rp) :: ur(lx, lx, lx)
 
 2075    real(kind=rp) :: us(lx, lx, lx)
 
 2076    real(kind=rp) :: ut(lx, lx, lx)
 
 2081          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 2082                    + dx(i,2) * u(2,j,1) &
 
 2083                    + dx(i,3) * u(3,j,1) &
 
 2084                    + dx(i,4) * u(4,j,1) &
 
 2085                    + dx(i,5) * u(5,j,1) &
 
 2086                    + dx(i,6) * u(6,j,1) &
 
 2087                    + dx(i,7) * u(7,j,1) &
 
 2088                    + dx(i,8) * u(8,j,1) &
 
 2089                    + dx(i,9) * u(9,j,1) &
 
 2090                    + dx(i,10) * u(10,j,1) &
 
 2091                    + dx(i,11) * u(11,j,1) &
 
 2092                    + dx(i,12) * u(12,j,1) &
 
 2093                    + dx(i,13) * u(13,j,1) &
 
 2094                    + dx(i,14) * u(14,j,1) &
 
 2095                    + dx(i,15) * u(15,j,1) &
 
 2096                    + dx(i,16) * u(16,j,1)
 
 2103             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2104                       + dy(j,2) * u(i,2,k) &
 
 2105                       + dy(j,3) * u(i,3,k) &
 
 2106                       + dy(j,4) * u(i,4,k) &
 
 2107                       + dy(j,5) * u(i,5,k) &
 
 2108                       + dy(j,6) * u(i,6,k) &
 
 2109                       + dy(j,7) * u(i,7,k) &
 
 2110                       + dy(j,8) * u(i,8,k) &
 
 2111                       + dy(j,9) * u(i,9,k) &
 
 2112                       + dy(j,10) * u(i,10,k) &
 
 2113                       + dy(j,11) * u(i,11,k) &
 
 2114                       + dy(j,12) * u(i,12,k) &
 
 2115                       + dy(j,13) * u(i,13,k) &
 
 2116                       + dy(j,14) * u(i,14,k) &
 
 2117                       + dy(j,15) * u(i,15,k) &
 
 2118                       + dy(j,16) * u(i,16,k)
 
 2125          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2126                    + dz(k,2) * u(i,1,2) &
 
 2127                    + dz(k,3) * u(i,1,3) &
 
 2128                    + dz(k,4) * u(i,1,4) &
 
 2129                    + dz(k,5) * u(i,1,5) &
 
 2130                    + dz(k,6) * u(i,1,6) &
 
 2131                    + dz(k,7) * u(i,1,7) &
 
 2132                    + dz(k,8) * u(i,1,8) &
 
 2133                    + dz(k,9) * u(i,1,9) &
 
 2134                    + dz(k,10) * u(i,1,10) &
 
 2135                    + dz(k,11) * u(i,1,11) &
 
 2136                    + dz(k,12) * u(i,1,12) &
 
 2137                    + dz(k,13) * u(i,1,13) &
 
 2138                    + dz(k,14) * u(i,1,14) &
 
 2139                    + dz(k,15) * u(i,1,15) &
 
 2140                    + dz(k,16) * u(i,1,16)
 
 2144    do i = 1, lx * lx * lx
 
 2145       ux(i,1,1) = w3(i,1,1) &
 
 2146                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 2147                     + dsdx(i,1,1) * us(i,1,1) &
 
 2148                     + dtdx(i,1,1) * ut(i,1,1) )
 
 2149       uy(i,1,1) = w3(i,1,1) &
 
 2150                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 2151                     + drdy(i,1,1) * ur(i,1,1) &
 
 2152                     + dtdy(i,1,1) * ut(i,1,1) )
 
 2153       uz(i,1,1) = w3(i,1,1) &
 
 2154                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 2155                     + drdz(i,1,1) * ur(i,1,1) &
 
 2156                     + dsdz(i,1,1) * us(i,1,1) )
 
 2158  end subroutine cpu_opgrad_lx16_single
 
 2160  subroutine cpu_opgrad_lx15_single(ux, uy, uz, u, dx, dy, dz, &
 
 2161       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 2162    integer, 
parameter :: lx = 15
 
 2163    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 2164    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 2165    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 2166    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 2167    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 2168    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 2169    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 2170    real(kind=rp) :: ur(lx, lx, lx)
 
 2171    real(kind=rp) :: us(lx, lx, lx)
 
 2172    real(kind=rp) :: ut(lx, lx, lx)
 
 2177          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 2178                    + dx(i,2) * u(2,j,1) &
 
 2179                    + dx(i,3) * u(3,j,1) &
 
 2180                    + dx(i,4) * u(4,j,1) &
 
 2181                    + dx(i,5) * u(5,j,1) &
 
 2182                    + dx(i,6) * u(6,j,1) &
 
 2183                    + dx(i,7) * u(7,j,1) &
 
 2184                    + dx(i,8) * u(8,j,1) &
 
 2185                    + dx(i,9) * u(9,j,1) &
 
 2186                    + dx(i,10) * u(10,j,1) &
 
 2187                    + dx(i,11) * u(11,j,1) &
 
 2188                    + dx(i,12) * u(12,j,1) &
 
 2189                    + dx(i,13) * u(13,j,1) &
 
 2190                    + dx(i,14) * u(14,j,1) &
 
 2191                    + dx(i,15) * u(15,j,1)
 
 2198             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2199                       + dy(j,2) * u(i,2,k) &
 
 2200                       + dy(j,3) * u(i,3,k) &
 
 2201                       + dy(j,4) * u(i,4,k) &
 
 2202                       + dy(j,5) * u(i,5,k) &
 
 2203                       + dy(j,6) * u(i,6,k) &
 
 2204                       + dy(j,7) * u(i,7,k) &
 
 2205                       + dy(j,8) * u(i,8,k) &
 
 2206                       + dy(j,9) * u(i,9,k) &
 
 2207                       + dy(j,10) * u(i,10,k) &
 
 2208                       + dy(j,11) * u(i,11,k) &
 
 2209                       + dy(j,12) * u(i,12,k) &
 
 2210                       + dy(j,13) * u(i,13,k) &
 
 2211                       + dy(j,14) * u(i,14,k) &
 
 2212                       + dy(j,15) * u(i,15,k)
 
 2219          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2220                    + dz(k,2) * u(i,1,2) &
 
 2221                    + dz(k,3) * u(i,1,3) &
 
 2222                    + dz(k,4) * u(i,1,4) &
 
 2223                    + dz(k,5) * u(i,1,5) &
 
 2224                    + dz(k,6) * u(i,1,6) &
 
 2225                    + dz(k,7) * u(i,1,7) &
 
 2226                    + dz(k,8) * u(i,1,8) &
 
 2227                    + dz(k,9) * u(i,1,9) &
 
 2228                    + dz(k,10) * u(i,1,10) &
 
 2229                    + dz(k,11) * u(i,1,11) &
 
 2230                    + dz(k,12) * u(i,1,12) &
 
 2231                    + dz(k,13) * u(i,1,13) &
 
 2232                    + dz(k,14) * u(i,1,14) &
 
 2233                    + dz(k,15) * u(i,1,15)
 
 2237    do i = 1, lx * lx * lx
 
 2238       ux(i,1,1) = w3(i,1,1) &
 
 2239                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 2240                     + dsdx(i,1,1) * us(i,1,1) &
 
 2241                     + dtdx(i,1,1) * ut(i,1,1) )
 
 2242       uy(i,1,1) = w3(i,1,1) &
 
 2243                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 2244                     + drdy(i,1,1) * ur(i,1,1) &
 
 2245                     + dtdy(i,1,1) * ut(i,1,1) )
 
 2246       uz(i,1,1) = w3(i,1,1) &
 
 2247                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 2248                     + drdz(i,1,1) * ur(i,1,1) &
 
 2249                     + dsdz(i,1,1) * us(i,1,1) )
 
 2251  end subroutine cpu_opgrad_lx15_single
 
 2253  subroutine cpu_opgrad_lx14_single(ux, uy, uz, u, dx, dy, dz, &
 
 2254       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 2255    integer, 
parameter :: lx = 14
 
 2256    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 2257    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 2258    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 2259    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 2260    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 2261    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 2262    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 2263    real(kind=rp) :: ur(lx, lx, lx)
 
 2264    real(kind=rp) :: us(lx, lx, lx)
 
 2265    real(kind=rp) :: ut(lx, lx, lx)
 
 2270          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 2271                    + dx(i,2) * u(2,j,1) &
 
 2272                    + dx(i,3) * u(3,j,1) &
 
 2273                    + dx(i,4) * u(4,j,1) &
 
 2274                    + dx(i,5) * u(5,j,1) &
 
 2275                    + dx(i,6) * u(6,j,1) &
 
 2276                    + dx(i,7) * u(7,j,1) &
 
 2277                    + dx(i,8) * u(8,j,1) &
 
 2278                    + dx(i,9) * u(9,j,1) &
 
 2279                    + dx(i,10) * u(10,j,1) &
 
 2280                    + dx(i,11) * u(11,j,1) &
 
 2281                    + dx(i,12) * u(12,j,1) &
 
 2282                    + dx(i,13) * u(13,j,1) &
 
 2283                    + dx(i,14) * u(14,j,1)
 
 2290             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2291                       + dy(j,2) * u(i,2,k) &
 
 2292                       + dy(j,3) * u(i,3,k) &
 
 2293                       + dy(j,4) * u(i,4,k) &
 
 2294                       + dy(j,5) * u(i,5,k) &
 
 2295                       + dy(j,6) * u(i,6,k) &
 
 2296                       + dy(j,7) * u(i,7,k) &
 
 2297                       + dy(j,8) * u(i,8,k) &
 
 2298                       + dy(j,9) * u(i,9,k) &
 
 2299                       + dy(j,10) * u(i,10,k) &
 
 2300                       + dy(j,11) * u(i,11,k) &
 
 2301                       + dy(j,12) * u(i,12,k) &
 
 2302                       + dy(j,13) * u(i,13,k) &
 
 2303                       + dy(j,14) * u(i,14,k)
 
 2310          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2311                    + dz(k,2) * u(i,1,2) &
 
 2312                    + dz(k,3) * u(i,1,3) &
 
 2313                    + dz(k,4) * u(i,1,4) &
 
 2314                    + dz(k,5) * u(i,1,5) &
 
 2315                    + dz(k,6) * u(i,1,6) &
 
 2316                    + dz(k,7) * u(i,1,7) &
 
 2317                    + dz(k,8) * u(i,1,8) &
 
 2318                    + dz(k,9) * u(i,1,9) &
 
 2319                    + dz(k,10) * u(i,1,10) &
 
 2320                    + dz(k,11) * u(i,1,11) &
 
 2321                    + dz(k,12) * u(i,1,12) &
 
 2322                    + dz(k,13) * u(i,1,13) &
 
 2323                    + dz(k,14) * u(i,1,14)
 
 2327    do i = 1, lx * lx * lx
 
 2328       ux(i,1,1) = w3(i,1,1) &
 
 2329                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 2330                     + dsdx(i,1,1) * us(i,1,1) &
 
 2331                     + dtdx(i,1,1) * ut(i,1,1) )
 
 2332       uy(i,1,1) = w3(i,1,1) &
 
 2333                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 2334                     + drdy(i,1,1) * ur(i,1,1) &
 
 2335                     + dtdy(i,1,1) * ut(i,1,1) )
 
 2336       uz(i,1,1) = w3(i,1,1) &
 
 2337                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 2338                     + drdz(i,1,1) * ur(i,1,1) &
 
 2339                     + dsdz(i,1,1) * us(i,1,1) )
 
 2341  end subroutine cpu_opgrad_lx14_single
 
 2343  subroutine cpu_opgrad_lx13_single(ux, uy, uz, u, dx, dy, dz, &
 
 2344       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 2345    integer, 
parameter :: lx = 13
 
 2346      real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 2347    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 2348    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 2349    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 2350    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 2351    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 2352    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 2353    real(kind=rp) :: ur(lx, lx, lx)
 
 2354    real(kind=rp) :: us(lx, lx, lx)
 
 2355    real(kind=rp) :: ut(lx, lx, lx)
 
 2360          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 2361                    + dx(i,2) * u(2,j,1) &
 
 2362                    + dx(i,3) * u(3,j,1) &
 
 2363                    + dx(i,4) * u(4,j,1) &
 
 2364                    + dx(i,5) * u(5,j,1) &
 
 2365                    + dx(i,6) * u(6,j,1) &
 
 2366                    + dx(i,7) * u(7,j,1) &
 
 2367                    + dx(i,8) * u(8,j,1) &
 
 2368                    + dx(i,9) * u(9,j,1) &
 
 2369                    + dx(i,10) * u(10,j,1) &
 
 2370                    + dx(i,11) * u(11,j,1) &
 
 2371                    + dx(i,12) * u(12,j,1) &
 
 2372                    + dx(i,13) * u(13,j,1)
 
 2379             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2380                       + dy(j,2) * u(i,2,k) &
 
 2381                       + dy(j,3) * u(i,3,k) &
 
 2382                       + dy(j,4) * u(i,4,k) &
 
 2383                       + dy(j,5) * u(i,5,k) &
 
 2384                       + dy(j,6) * u(i,6,k) &
 
 2385                       + dy(j,7) * u(i,7,k) &
 
 2386                       + dy(j,8) * u(i,8,k) &
 
 2387                       + dy(j,9) * u(i,9,k) &
 
 2388                       + dy(j,10) * u(i,10,k) &
 
 2389                       + dy(j,11) * u(i,11,k) &
 
 2390                       + dy(j,12) * u(i,12,k) &
 
 2391                       + dy(j,13) * u(i,13,k)
 
 2398          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2399                    + dz(k,2) * u(i,1,2) &
 
 2400                    + dz(k,3) * u(i,1,3) &
 
 2401                    + dz(k,4) * u(i,1,4) &
 
 2402                    + dz(k,5) * u(i,1,5) &
 
 2403                    + dz(k,6) * u(i,1,6) &
 
 2404                    + dz(k,7) * u(i,1,7) &
 
 2405                    + dz(k,8) * u(i,1,8) &
 
 2406                    + dz(k,9) * u(i,1,9) &
 
 2407                    + dz(k,10) * u(i,1,10) &
 
 2408                    + dz(k,11) * u(i,1,11) &
 
 2409                    + dz(k,12) * u(i,1,12) &
 
 2410                    + dz(k,13) * u(i,1,13)
 
 2414    do i = 1, lx * lx * lx
 
 2415       ux(i,1,1) = w3(i,1,1) &
 
 2416                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 2417                     + dsdx(i,1,1) * us(i,1,1) &
 
 2418                     + dtdx(i,1,1) * ut(i,1,1) )
 
 2419       uy(i,1,1) = w3(i,1,1) &
 
 2420                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 2421                     + drdy(i,1,1) * ur(i,1,1) &
 
 2422                     + dtdy(i,1,1) * ut(i,1,1) )
 
 2423       uz(i,1,1) = w3(i,1,1) &
 
 2424                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 2425                     + drdz(i,1,1) * ur(i,1,1) &
 
 2426                     + dsdz(i,1,1) * us(i,1,1) )
 
 2428  end subroutine cpu_opgrad_lx13_single
 
 2430  subroutine cpu_opgrad_lx12_single(ux, uy, uz, u, dx, dy, dz, &
 
 2431       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 2432    integer, 
parameter :: lx = 12
 
 2433    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 2434    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 2435    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 2436    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 2437    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 2438    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 2439    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 2440    real(kind=rp) :: ur(lx, lx, lx)
 
 2441    real(kind=rp) :: us(lx, lx, lx)
 
 2442    real(kind=rp) :: ut(lx, lx, lx)
 
 2447          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 2448                    + dx(i,2) * u(2,j,1) &
 
 2449                    + dx(i,3) * u(3,j,1) &
 
 2450                    + dx(i,4) * u(4,j,1) &
 
 2451                    + dx(i,5) * u(5,j,1) &
 
 2452                    + dx(i,6) * u(6,j,1) &
 
 2453                    + dx(i,7) * u(7,j,1) &
 
 2454                    + dx(i,8) * u(8,j,1) &
 
 2455                    + dx(i,9) * u(9,j,1) &
 
 2456                    + dx(i,10) * u(10,j,1) &
 
 2457                    + dx(i,11) * u(11,j,1) &
 
 2458                    + dx(i,12) * u(12,j,1)
 
 2465             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2466                       + dy(j,2) * u(i,2,k) &
 
 2467                       + dy(j,3) * u(i,3,k) &
 
 2468                       + dy(j,4) * u(i,4,k) &
 
 2469                       + dy(j,5) * u(i,5,k) &
 
 2470                       + dy(j,6) * u(i,6,k) &
 
 2471                       + dy(j,7) * u(i,7,k) &
 
 2472                       + dy(j,8) * u(i,8,k) &
 
 2473                       + dy(j,9) * u(i,9,k) &
 
 2474                       + dy(j,10) * u(i,10,k) &
 
 2475                       + dy(j,11) * u(i,11,k) &
 
 2476                       + dy(j,12) * u(i,12,k)
 
 2483          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2484                    + dz(k,2) * u(i,1,2) &
 
 2485                    + dz(k,3) * u(i,1,3) &
 
 2486                    + dz(k,4) * u(i,1,4) &
 
 2487                    + dz(k,5) * u(i,1,5) &
 
 2488                    + dz(k,6) * u(i,1,6) &
 
 2489                    + dz(k,7) * u(i,1,7) &
 
 2490                    + dz(k,8) * u(i,1,8) &
 
 2491                    + dz(k,9) * u(i,1,9) &
 
 2492                    + dz(k,10) * u(i,1,10) &
 
 2493                    + dz(k,11) * u(i,1,11) &
 
 2494                    + dz(k,12) * u(i,1,12)
 
 2498    do i = 1, lx * lx * lx
 
 2499       ux(i,1,1) = w3(i,1,1) &
 
 2500                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 2501                     + dsdx(i,1,1) * us(i,1,1) &
 
 2502                     + dtdx(i,1,1) * ut(i,1,1) )
 
 2503       uy(i,1,1) = w3(i,1,1) &
 
 2504                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 2505                     + drdy(i,1,1) * ur(i,1,1) &
 
 2506                     + dtdy(i,1,1) * ut(i,1,1) )
 
 2507       uz(i,1,1) = w3(i,1,1) &
 
 2508                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 2509                     + drdz(i,1,1) * ur(i,1,1) &
 
 2510                     + dsdz(i,1,1) * us(i,1,1) )
 
 2512  end subroutine cpu_opgrad_lx12_single
 
 2514  subroutine cpu_opgrad_lx11_single(ux, uy, uz, u, dx, dy, dz, &
 
 2515       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 2516    integer, 
parameter :: lx = 11
 
 2517    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 2518    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 2519    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 2520    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 2521    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 2522    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 2523    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 2524    real(kind=rp) :: ur(lx, lx, lx)
 
 2525    real(kind=rp) :: us(lx, lx, lx)
 
 2526    real(kind=rp) :: ut(lx, lx, lx)
 
 2531          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 2532                    + dx(i,2) * u(2,j,1) &
 
 2533                    + dx(i,3) * u(3,j,1) &
 
 2534                    + dx(i,4) * u(4,j,1) &
 
 2535                    + dx(i,5) * u(5,j,1) &
 
 2536                    + dx(i,6) * u(6,j,1) &
 
 2537                    + dx(i,7) * u(7,j,1) &
 
 2538                    + dx(i,8) * u(8,j,1) &
 
 2539                    + dx(i,9) * u(9,j,1) &
 
 2540                    + dx(i,10) * u(10,j,1) &
 
 2541                    + dx(i,11) * u(11,j,1)
 
 2548             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2549                       + dy(j,2) * u(i,2,k) &
 
 2550                       + dy(j,3) * u(i,3,k) &
 
 2551                       + dy(j,4) * u(i,4,k) &
 
 2552                       + dy(j,5) * u(i,5,k) &
 
 2553                       + dy(j,6) * u(i,6,k) &
 
 2554                       + dy(j,7) * u(i,7,k) &
 
 2555                       + dy(j,8) * u(i,8,k) &
 
 2556                       + dy(j,9) * u(i,9,k) &
 
 2557                       + dy(j,10) * u(i,10,k) &
 
 2558                       + dy(j,11) * u(i,11,k)
 
 2565          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2566                    + dz(k,2) * u(i,1,2) &
 
 2567                    + dz(k,3) * u(i,1,3) &
 
 2568                    + dz(k,4) * u(i,1,4) &
 
 2569                    + dz(k,5) * u(i,1,5) &
 
 2570                    + dz(k,6) * u(i,1,6) &
 
 2571                    + dz(k,7) * u(i,1,7) &
 
 2572                    + dz(k,8) * u(i,1,8) &
 
 2573                    + dz(k,9) * u(i,1,9) &
 
 2574                    + dz(k,10) * u(i,1,10) &
 
 2575                    + dz(k,11) * u(i,1,11)
 
 2579    do i = 1, lx * lx * lx
 
 2580       ux(i,1,1) = w3(i,1,1) &
 
 2581                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 2582                     + dsdx(i,1,1) * us(i,1,1) &
 
 2583                     + dtdx(i,1,1) * ut(i,1,1) )
 
 2584       uy(i,1,1) = w3(i,1,1) &
 
 2585                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 2586                     + drdy(i,1,1) * ur(i,1,1) &
 
 2587                     + dtdy(i,1,1) * ut(i,1,1) )
 
 2588       uz(i,1,1) = w3(i,1,1) &
 
 2589                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 2590                     + drdz(i,1,1) * ur(i,1,1) &
 
 2591                     + dsdz(i,1,1) * us(i,1,1) )
 
 2593  end subroutine cpu_opgrad_lx11_single
 
 2595  subroutine cpu_opgrad_lx10_single(ux, uy, uz, u, dx, dy, dz, &
 
 2596       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 2597    integer, 
parameter :: lx = 10
 
 2598    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 2599    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 2600    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 2601    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 2602    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 2603    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 2604    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 2605    real(kind=rp) :: ur(lx, lx, lx)
 
 2606    real(kind=rp) :: us(lx, lx, lx)
 
 2607    real(kind=rp) :: ut(lx, lx, lx)
 
 2612          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 2613                    + dx(i,2) * u(2,j,1) &
 
 2614                    + dx(i,3) * u(3,j,1) &
 
 2615                    + dx(i,4) * u(4,j,1) &
 
 2616                    + dx(i,5) * u(5,j,1) &
 
 2617                    + dx(i,6) * u(6,j,1) &
 
 2618                    + dx(i,7) * u(7,j,1) &
 
 2619                    + dx(i,8) * u(8,j,1) &
 
 2620                    + dx(i,9) * u(9,j,1) &
 
 2621                    + dx(i,10) * u(10,j,1)
 
 2628             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2629                       + dy(j,2) * u(i,2,k) &
 
 2630                       + dy(j,3) * u(i,3,k) &
 
 2631                       + dy(j,4) * u(i,4,k) &
 
 2632                       + dy(j,5) * u(i,5,k) &
 
 2633                       + dy(j,6) * u(i,6,k) &
 
 2634                       + dy(j,7) * u(i,7,k) &
 
 2635                       + dy(j,8) * u(i,8,k) &
 
 2636                       + dy(j,9) * u(i,9,k) &
 
 2637                       + dy(j,10) * u(i,10,k)
 
 2644          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2645                    + dz(k,2) * u(i,1,2) &
 
 2646                    + dz(k,3) * u(i,1,3) &
 
 2647                    + dz(k,4) * u(i,1,4) &
 
 2648                    + dz(k,5) * u(i,1,5) &
 
 2649                    + dz(k,6) * u(i,1,6) &
 
 2650                    + dz(k,7) * u(i,1,7) &
 
 2651                    + dz(k,8) * u(i,1,8) &
 
 2652                    + dz(k,9) * u(i,1,9) &
 
 2653                    + dz(k,10) * u(i,1,10)
 
 2657    do i = 1, lx * lx * lx
 
 2658       ux(i,1,1) = w3(i,1,1) &
 
 2659                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 2660                     + dsdx(i,1,1) * us(i,1,1) &
 
 2661                     + dtdx(i,1,1) * ut(i,1,1) )
 
 2662       uy(i,1,1) = w3(i,1,1) &
 
 2663                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 2664                     + drdy(i,1,1) * ur(i,1,1) &
 
 2665                     + dtdy(i,1,1) * ut(i,1,1) )
 
 2666       uz(i,1,1) = w3(i,1,1) &
 
 2667                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 2668                     + drdz(i,1,1) * ur(i,1,1) &
 
 2669                     + dsdz(i,1,1) * us(i,1,1) )
 
 2671  end subroutine cpu_opgrad_lx10_single
 
 2673  subroutine cpu_opgrad_lx9_single(ux, uy, uz, u, dx, dy, dz, &
 
 2674       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 2675    integer, 
parameter :: lx = 9
 
 2676    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 2677    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 2678    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 2679    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 2680    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 2681    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 2682    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 2683    real(kind=rp) :: ur(lx, lx, lx)
 
 2684    real(kind=rp) :: us(lx, lx, lx)
 
 2685    real(kind=rp) :: ut(lx, lx, lx)
 
 2690          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 2691                    + dx(i,2) * u(2,j,1) &
 
 2692                    + dx(i,3) * u(3,j,1) &
 
 2693                    + dx(i,4) * u(4,j,1) &
 
 2694                    + dx(i,5) * u(5,j,1) &
 
 2695                    + dx(i,6) * u(6,j,1) &
 
 2696                    + dx(i,7) * u(7,j,1) &
 
 2697                    + dx(i,8) * u(8,j,1) &
 
 2698                    + dx(i,9) * u(9,j,1)
 
 2705             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2706                       + dy(j,2) * u(i,2,k) &
 
 2707                       + dy(j,3) * u(i,3,k) &
 
 2708                       + dy(j,4) * u(i,4,k) &
 
 2709                       + dy(j,5) * u(i,5,k) &
 
 2710                       + dy(j,6) * u(i,6,k) &
 
 2711                       + dy(j,7) * u(i,7,k) &
 
 2712                       + dy(j,8) * u(i,8,k) &
 
 2713                       + dy(j,9) * u(i,9,k)
 
 2720          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2721                    + dz(k,2) * u(i,1,2) &
 
 2722                    + dz(k,3) * u(i,1,3) &
 
 2723                    + dz(k,4) * u(i,1,4) &
 
 2724                    + dz(k,5) * u(i,1,5) &
 
 2725                    + dz(k,6) * u(i,1,6) &
 
 2726                    + dz(k,7) * u(i,1,7) &
 
 2727                    + dz(k,8) * u(i,1,8) &
 
 2728                    + dz(k,9) * u(i,1,9)
 
 2732    do i = 1, lx * lx * lx
 
 2733       ux(i,1,1) = w3(i,1,1) &
 
 2734                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 2735                     + dsdx(i,1,1) * us(i,1,1) &
 
 2736                     + dtdx(i,1,1) * ut(i,1,1) )
 
 2737       uy(i,1,1) = w3(i,1,1) &
 
 2738                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 2739                     + drdy(i,1,1) * ur(i,1,1) &
 
 2740                     + dtdy(i,1,1) * ut(i,1,1) )
 
 2741       uz(i,1,1) = w3(i,1,1) &
 
 2742                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 2743                     + drdz(i,1,1) * ur(i,1,1) &
 
 2744                     + dsdz(i,1,1) * us(i,1,1) )
 
 2746  end subroutine cpu_opgrad_lx9_single
 
 2748  subroutine cpu_opgrad_lx8_single(ux, uy, uz, u, dx, dy, dz, &
 
 2749       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 2750    integer, 
parameter :: lx = 8
 
 2751    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 2752    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 2753    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 2754    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 2755    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 2756    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 2757    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 2758    real(kind=rp) :: ur(lx, lx, lx)
 
 2759    real(kind=rp) :: us(lx, lx, lx)
 
 2760    real(kind=rp) :: ut(lx, lx, lx)
 
 2765          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 2766                    + dx(i,2) * u(2,j,1) &
 
 2767                    + dx(i,3) * u(3,j,1) &
 
 2768                    + dx(i,4) * u(4,j,1) &
 
 2769                    + dx(i,5) * u(5,j,1) &
 
 2770                    + dx(i,6) * u(6,j,1) &
 
 2771                    + dx(i,7) * u(7,j,1) &
 
 2772                    + dx(i,8) * u(8,j,1)
 
 2779             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2780                       + dy(j,2) * u(i,2,k) &
 
 2781                       + dy(j,3) * u(i,3,k) &
 
 2782                       + dy(j,4) * u(i,4,k) &
 
 2783                       + dy(j,5) * u(i,5,k) &
 
 2784                       + dy(j,6) * u(i,6,k) &
 
 2785                       + dy(j,7) * u(i,7,k) &
 
 2786                       + dy(j,8) * u(i,8,k)
 
 2793          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2794                    + dz(k,2) * u(i,1,2) &
 
 2795                    + dz(k,3) * u(i,1,3) &
 
 2796                    + dz(k,4) * u(i,1,4) &
 
 2797                    + dz(k,5) * u(i,1,5) &
 
 2798                    + dz(k,6) * u(i,1,6) &
 
 2799                    + dz(k,7) * u(i,1,7) &
 
 2800                    + dz(k,8) * u(i,1,8)
 
 2804    do i = 1, lx * lx * lx
 
 2805       ux(i,1,1) = w3(i,1,1) &
 
 2806                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 2807                     + dsdx(i,1,1) * us(i,1,1) &
 
 2808                     + dtdx(i,1,1) * ut(i,1,1) )
 
 2809       uy(i,1,1) = w3(i,1,1) &
 
 2810                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 2811                     + drdy(i,1,1) * ur(i,1,1) &
 
 2812                     + dtdy(i,1,1) * ut(i,1,1) )
 
 2813       uz(i,1,1) = w3(i,1,1) &
 
 2814                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 2815                     + drdz(i,1,1) * ur(i,1,1) &
 
 2816                     + dsdz(i,1,1) * us(i,1,1) )
 
 2818  end subroutine cpu_opgrad_lx8_single
 
 2820  subroutine cpu_opgrad_lx7_single(ux, uy, uz, u, dx, dy, dz, &
 
 2821       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 2822    integer, 
parameter :: lx = 7
 
 2823    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 2824    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 2825    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 2826    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 2827    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 2828    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 2829    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 2830    real(kind=rp) :: ur(lx, lx, lx)
 
 2831    real(kind=rp) :: us(lx, lx, lx)
 
 2832    real(kind=rp) :: ut(lx, lx, lx)
 
 2837          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 2838                    + dx(i,2) * u(2,j,1) &
 
 2839                    + dx(i,3) * u(3,j,1) &
 
 2840                    + dx(i,4) * u(4,j,1) &
 
 2841                    + dx(i,5) * u(5,j,1) &
 
 2842                    + dx(i,6) * u(6,j,1) &
 
 2843                    + dx(i,7) * u(7,j,1)
 
 2850             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2851                       + dy(j,2) * u(i,2,k) &
 
 2852                       + dy(j,3) * u(i,3,k) &
 
 2853                       + dy(j,4) * u(i,4,k) &
 
 2854                       + dy(j,5) * u(i,5,k) &
 
 2855                       + dy(j,6) * u(i,6,k) &
 
 2856                       + dy(j,7) * u(i,7,k)
 
 2863          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2864                    + dz(k,2) * u(i,1,2) &
 
 2865                    + dz(k,3) * u(i,1,3) &
 
 2866                    + dz(k,4) * u(i,1,4) &
 
 2867                    + dz(k,5) * u(i,1,5) &
 
 2868                    + dz(k,6) * u(i,1,6) &
 
 2869                    + dz(k,7) * u(i,1,7)
 
 2873    do i = 1, lx * lx * lx
 
 2874       ux(i,1,1) = w3(i,1,1) &
 
 2875                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 2876                     + dsdx(i,1,1) * us(i,1,1) &
 
 2877                     + dtdx(i,1,1) * ut(i,1,1) )
 
 2878       uy(i,1,1) = w3(i,1,1) &
 
 2879                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 2880                     + drdy(i,1,1) * ur(i,1,1) &
 
 2881                     + dtdy(i,1,1) * ut(i,1,1) )
 
 2882       uz(i,1,1) = w3(i,1,1) &
 
 2883                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 2884                     + drdz(i,1,1) * ur(i,1,1) &
 
 2885                     + dsdz(i,1,1) * us(i,1,1) )
 
 2887  end subroutine cpu_opgrad_lx7_single
 
 2889  subroutine cpu_opgrad_lx6_single(ux, uy, uz, u, dx, dy, dz, &
 
 2890       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 2891    integer, 
parameter :: lx = 6
 
 2892    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 2893    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 2894    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 2895    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 2896    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 2897    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 2898    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 2899    real(kind=rp) :: ur(lx, lx, lx)
 
 2900    real(kind=rp) :: us(lx, lx, lx)
 
 2901    real(kind=rp) :: ut(lx, lx, lx)
 
 2906          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 2907                    + dx(i,2) * u(2,j,1) &
 
 2908                    + dx(i,3) * u(3,j,1) &
 
 2909                    + dx(i,4) * u(4,j,1) &
 
 2910                    + dx(i,5) * u(5,j,1) &
 
 2911                    + dx(i,6) * u(6,j,1)
 
 2918             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2919                       + dy(j,2) * u(i,2,k) &
 
 2920                       + dy(j,3) * u(i,3,k) &
 
 2921                       + dy(j,4) * u(i,4,k) &
 
 2922                       + dy(j,5) * u(i,5,k) &
 
 2923                       + dy(j,6) * u(i,6,k)
 
 2930          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2931                    + dz(k,2) * u(i,1,2) &
 
 2932                    + dz(k,3) * u(i,1,3) &
 
 2933                    + dz(k,4) * u(i,1,4) &
 
 2934                    + dz(k,5) * u(i,1,5) &
 
 2935                    + dz(k,6) * u(i,1,6)
 
 2939    do i = 1, lx * lx * lx
 
 2940       ux(i,1,1) = w3(i,1,1) &
 
 2941                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 2942                     + dsdx(i,1,1) * us(i,1,1) &
 
 2943                     + dtdx(i,1,1) * ut(i,1,1) )
 
 2944       uy(i,1,1) = w3(i,1,1) &
 
 2945                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 2946                     + drdy(i,1,1) * ur(i,1,1) &
 
 2947                     + dtdy(i,1,1) * ut(i,1,1) )
 
 2948       uz(i,1,1) = w3(i,1,1) &
 
 2949                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 2950                     + drdz(i,1,1) * ur(i,1,1) &
 
 2951                     + dsdz(i,1,1) * us(i,1,1) )
 
 2953  end subroutine cpu_opgrad_lx6_single
 
 2955  subroutine cpu_opgrad_lx5_single(ux, uy, uz, u, dx, dy, dz, &
 
 2956       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 2957    integer, 
parameter :: lx = 5
 
 2958    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 2959    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 2960    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 2961    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 2962    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 2963    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 2964    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 2965    real(kind=rp) :: ur(lx, lx, lx)
 
 2966    real(kind=rp) :: us(lx, lx, lx)
 
 2967    real(kind=rp) :: ut(lx, lx, lx)
 
 2972          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 2973                    + dx(i,2) * u(2,j,1) &
 
 2974                    + dx(i,3) * u(3,j,1) &
 
 2975                    + dx(i,4) * u(4,j,1) &
 
 2976                    + dx(i,5) * u(5,j,1)
 
 2983             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 2984                       + dy(j,2) * u(i,2,k) &
 
 2985                       + dy(j,3) * u(i,3,k) &
 
 2986                       + dy(j,4) * u(i,4,k) &
 
 2987                       + dy(j,5) * u(i,5,k)
 
 2994          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 2995                    + dz(k,2) * u(i,1,2) &
 
 2996                    + dz(k,3) * u(i,1,3) &
 
 2997                    + dz(k,4) * u(i,1,4) &
 
 2998                    + dz(k,5) * u(i,1,5)
 
 3002    do i = 1, lx * lx * lx
 
 3003       ux(i,1,1) = w3(i,1,1) &
 
 3004                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 3005                     + dsdx(i,1,1) * us(i,1,1) &
 
 3006                     + dtdx(i,1,1) * ut(i,1,1) )
 
 3007       uy(i,1,1) = w3(i,1,1) &
 
 3008                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 3009                     + drdy(i,1,1) * ur(i,1,1) &
 
 3010                     + dtdy(i,1,1) * ut(i,1,1) )
 
 3011       uz(i,1,1) = w3(i,1,1) &
 
 3012                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 3013                     + drdz(i,1,1) * ur(i,1,1) &
 
 3014                     + dsdz(i,1,1) * us(i,1,1) )
 
 3016  end subroutine cpu_opgrad_lx5_single
 
 3018  subroutine cpu_opgrad_lx4_single(ux, uy, uz, u, dx, dy, dz, &
 
 3019       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 3020    integer, 
parameter :: lx = 4
 
 3021    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 3022    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 3023    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 3024    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 3025    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 3026    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 3027    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 3028    real(kind=rp) :: ur(lx, lx, lx)
 
 3029    real(kind=rp) :: us(lx, lx, lx)
 
 3030    real(kind=rp) :: ut(lx, lx, lx)
 
 3035          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 3036                    + dx(i,2) * u(2,j,1) &
 
 3037                    + dx(i,3) * u(3,j,1) &
 
 3038                    + dx(i,4) * u(4,j,1)
 
 3045             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 3046                       + dy(j,2) * u(i,2,k) &
 
 3047                       + dy(j,3) * u(i,3,k) &
 
 3048                       + dy(j,4) * u(i,4,k)
 
 3055          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 3056                    + dz(k,2) * u(i,1,2) &
 
 3057                    + dz(k,3) * u(i,1,3) &
 
 3058                    + dz(k,4) * u(i,1,4)
 
 3062    do i = 1, lx * lx * lx
 
 3063       ux(i,1,1) = w3(i,1,1) &
 
 3064                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 3065                     + dsdx(i,1,1) * us(i,1,1) &
 
 3066                     + dtdx(i,1,1) * ut(i,1,1) )
 
 3067       uy(i,1,1) = w3(i,1,1) &
 
 3068                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 3069                     + drdy(i,1,1) * ur(i,1,1) &
 
 3070                     + dtdy(i,1,1) * ut(i,1,1) )
 
 3071       uz(i,1,1) = w3(i,1,1) &
 
 3072                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 3073                     + drdz(i,1,1) * ur(i,1,1) &
 
 3074                     + dsdz(i,1,1) * us(i,1,1) )
 
 3076  end subroutine cpu_opgrad_lx4_single
 
 3078  subroutine cpu_opgrad_lx3_single(ux, uy, uz, u, dx, dy, dz, &
 
 3079       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 3080    integer, 
parameter :: lx = 3
 
 3081    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 3082    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 3083    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 3084    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 3085    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 3086    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 3087    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 3088    real(kind=rp) :: ur(lx, lx, lx)
 
 3089    real(kind=rp) :: us(lx, lx, lx)
 
 3090    real(kind=rp) :: ut(lx, lx, lx)
 
 3095          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 3096                    + dx(i,2) * u(2,j,1) &
 
 3097                    + dx(i,3) * u(3,j,1)
 
 3104             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 3105                       + dy(j,2) * u(i,2,k) &
 
 3106                       + dy(j,3) * u(i,3,k)
 
 3113          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 3114                    + dz(k,2) * u(i,1,2) &
 
 3115                    + dz(k,3) * u(i,1,3)
 
 3119    do i = 1, lx * lx * lx
 
 3120       ux(i,1,1) = w3(i,1,1) &
 
 3121                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 3122                     + dsdx(i,1,1) * us(i,1,1) &
 
 3123                     + dtdx(i,1,1) * ut(i,1,1) )
 
 3124       uy(i,1,1) = w3(i,1,1) &
 
 3125                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 3126                     + drdy(i,1,1) * ur(i,1,1) &
 
 3127                     + dtdy(i,1,1) * ut(i,1,1) )
 
 3128       uz(i,1,1) = w3(i,1,1) &
 
 3129                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 3130                     + drdz(i,1,1) * ur(i,1,1) &
 
 3131                     + dsdz(i,1,1) * us(i,1,1) )
 
 3133  end subroutine cpu_opgrad_lx3_single
 
 3135  subroutine cpu_opgrad_lx2_single(ux, uy, uz, u, dx, dy, dz, &
 
 3136       drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
 
 3137    integer, 
parameter :: lx = 2
 
 3138    real(kind=rp), 
dimension(lx, lx, lx), 
intent(inout) :: ux, uy, uz
 
 3139    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: u
 
 3140    real(kind=rp), 
dimension(lx, lx), 
intent(in) :: dx, dy, dz
 
 3141    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdx, dsdx, dtdx
 
 3142    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdy, dsdy, dtdy
 
 3143    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: drdz, dsdz, dtdz
 
 3144    real(kind=rp), 
dimension(lx, lx, lx), 
intent(in) :: w3
 
 3145    real(kind=rp) :: ur(lx, lx, lx)
 
 3146    real(kind=rp) :: us(lx, lx, lx)
 
 3147    real(kind=rp) :: ut(lx, lx, lx)
 
 3152          ur(i,j,1) = dx(i,1) * u(1,j,1) &
 
 3153                    + dx(i,2) * u(2,j,1)
 
 3160             us(i,j,k) = dy(j,1) * u(i,1,k) &
 
 3161                       + dy(j,2) * u(i,2,k)
 
 3168          ut(i,1,k) = dz(k,1) * u(i,1,1) &
 
 3169                    + dz(k,2) * u(i,1,2)
 
 3173    do i = 1, lx * lx * lx
 
 3174       ux(i,1,1) = w3(i,1,1) &
 
 3175                   * ( drdx(i,1,1) * ur(i,1,1) &
 
 3176                     + dsdx(i,1,1) * us(i,1,1) &
 
 3177                     + dtdx(i,1,1) * ut(i,1,1) )
 
 3178       uy(i,1,1) = w3(i,1,1) &
 
 3179                   * ( dsdy(i,1,1) * us(i,1,1) &
 
 3180                     + drdy(i,1,1) * ur(i,1,1) &
 
 3181                     + dtdy(i,1,1) * ut(i,1,1) )
 
 3182       uz(i,1,1) = w3(i,1,1) &
 
 3183                   * ( dtdz(i,1,1) * ut(i,1,1) &
 
 3184                     + drdz(i,1,1) * ur(i,1,1) &
 
 3185                     + dsdz(i,1,1) * us(i,1,1) )
 
 3188  end subroutine cpu_opgrad_lx2_single
 
 3190end submodule cpu_opgrad