38    type(
field_t), 
intent(inout) :: u,v, w
 
   39    type(
field_t), 
intent(inout) :: uu, vv, ww
 
   41    real(kind=
rp), 
dimension(3), 
intent(in) :: ab
 
   42    integer, 
intent(in) :: nab
 
   47    do concurrent(i = 1:n)
 
   48       u%x(i,1,1,1) = ab(1) * uu%x(i,1,1,1) + ab(2) * uulag%lf(1)%x(i,1,1,1)
 
   49       v%x(i,1,1,1) = ab(1) * vv%x(i,1,1,1) + ab(2) * vvlag%lf(1)%x(i,1,1,1)
 
   50       w%x(i,1,1,1) = ab(1) * ww%x(i,1,1,1) + ab(2) * wwlag%lf(1)%x(i,1,1,1)
 
   54       do concurrent(i = 1:n)
 
   55          u%x(i,1,1,1) = u%x(i,1,1,1) + ab(3) * uulag%lf(2)%x(i,1,1,1)
 
   56          v%x(i,1,1,1) = v%x(i,1,1,1) + ab(3) * vvlag%lf(2)%x(i,1,1,1)
 
   57          w%x(i,1,1,1) = w%x(i,1,1,1) + ab(3) * wwlag%lf(2)%x(i,1,1,1)
 
 
   64                             fx_laglag, fy_laglag, fz_laglag, fx, fy, fz, &
 
   66    type(
field_t), 
intent(inout) :: fx_lag, fy_lag, fz_lag
 
   67    type(
field_t), 
intent(inout) :: fx_laglag, fy_laglag, fz_laglag
 
   68    real(kind=
rp), 
intent(inout) :: rho, ext_coeffs(4)
 
   69    integer, 
intent(in) :: n
 
   70    real(kind=
rp), 
intent(inout) :: fx(n), fy(n), fz(n)
 
   72    type(
field_t), 
pointer :: temp1, temp2, temp3
 
   73    integer :: temp_indices(3)
 
   79    do concurrent(i = 1:n)
 
   80       temp1%x(i,1,1,1) = ext_coeffs(2) * fx_lag%x(i,1,1,1) + &
 
   81                          ext_coeffs(3) * fx_laglag%x(i,1,1,1)
 
   82       temp2%x(i,1,1,1) = ext_coeffs(2) * fy_lag%x(i,1,1,1) + &
 
   83                          ext_coeffs(3) * fy_laglag%x(i,1,1,1)
 
   84       temp3%x(i,1,1,1) = ext_coeffs(2) * fz_lag%x(i,1,1,1) + &
 
   85                          ext_coeffs(3) * fz_laglag%x(i,1,1,1)
 
   88    do concurrent(i = 1:n)
 
   89       fx_laglag%x(i,1,1,1) = fx_lag%x(i,1,1,1)
 
   90       fy_laglag%x(i,1,1,1) = fy_lag%x(i,1,1,1)
 
   91       fz_laglag%x(i,1,1,1) = fz_lag%x(i,1,1,1)
 
   92       fx_lag%x(i,1,1,1) = fx(i)
 
   93       fy_lag%x(i,1,1,1) = fy(i)
 
   94       fz_lag%x(i,1,1,1) = fz(i)
 
   97    do concurrent(i = 1:n)
 
   98       fx(i) = (ext_coeffs(1) * fx(i) + temp1%x(i,1,1,1)) * rho
 
   99       fy(i) = (ext_coeffs(1) * fy(i) + temp2%x(i,1,1,1)) * rho
 
  100       fz(i) = (ext_coeffs(1) * fz(i) + temp3%x(i,1,1,1)) * rho
 
 
  109    type(
field_t), 
intent(inout) :: fs_lag
 
  110    type(
field_t), 
intent(inout) :: fs_laglag
 
  111    real(kind=
rp), 
intent(inout) :: rho, ext_coeffs(4)
 
  112    integer, 
intent(in) :: n
 
  113    real(kind=
rp), 
intent(inout) :: fs(n)
 
  115    type(
field_t), 
pointer :: temp1
 
  116    integer :: temp_index
 
  120    do concurrent(i = 1:n)
 
  121       temp1%x(i,1,1,1) = ext_coeffs(2) * fs_lag%x(i,1,1,1) + &
 
  122                          ext_coeffs(3) * fs_laglag%x(i,1,1,1)
 
  125    do concurrent(i = 1:n)
 
  126       fs_laglag%x(i,1,1,1) = fs_lag%x(i,1,1,1)
 
  127       fs_lag%x(i,1,1,1) = fs(i)
 
  130    do concurrent(i = 1:n)
 
  131       fs(i) = (ext_coeffs(1) * fs(i) + temp1%x(i,1,1,1)) * rho
 
 
  138                               u, v, w, B, rho, dt, bd, nbd, n)
 
  139    integer, 
intent(in) :: n, nbd
 
  140    type(
field_t), 
intent(in) :: u, v, w
 
  142    real(kind=
rp), 
intent(inout) :: bfx(n), bfy(n), bfz(n)
 
  143    real(kind=
rp), 
intent(in) :: b(n)
 
  144    real(kind=
rp), 
intent(in) :: dt, rho, bd(4)
 
  145    type(
field_t), 
pointer :: tb1, tb2, tb3
 
  146    type(
field_t), 
pointer :: ta1, ta2, ta3
 
  147    integer :: temp_indices(6)
 
  157    do concurrent(i = 1:n)
 
  158       tb1%x(i,1,1,1) = u%x(i,1,1,1) * b(i) * bd(2)
 
  159       tb2%x(i,1,1,1) = v%x(i,1,1,1) * b(i) * bd(2)
 
  160       tb3%x(i,1,1,1) = w%x(i,1,1,1) * b(i) * bd(2)
 
  164       do concurrent(i = 1:n)
 
  165          ta1%x(i,1,1,1) = ulag%lf(ilag-1)%x(i,1,1,1) * b(i) * bd(ilag+1)
 
  166          ta2%x(i,1,1,1) = vlag%lf(ilag-1)%x(i,1,1,1) * b(i) * bd(ilag+1)
 
  167          ta3%x(i,1,1,1) = wlag%lf(ilag-1)%x(i,1,1,1) * b(i) * bd(ilag+1)
 
  170       do concurrent(i = 1:n)
 
  171          tb1%x(i,1,1,1) = tb1%x(i,1,1,1) + ta1%x(i,1,1,1)
 
  172          tb2%x(i,1,1,1) = tb2%x(i,1,1,1) + ta2%x(i,1,1,1)
 
  173          tb3%x(i,1,1,1) = tb3%x(i,1,1,1) + ta3%x(i,1,1,1)
 
  177    do concurrent(i = 1:n)
 
  178       bfx(i) = bfx(i) + tb1%x(i,1,1,1) * (rho / dt)
 
  179       bfy(i) = bfy(i) + tb2%x(i,1,1,1) * (rho / dt)
 
  180       bfz(i) = bfz(i) + tb3%x(i,1,1,1) * (rho / dt)
 
 
  188    integer, 
intent(in) :: n, nbd
 
  191    real(kind=
rp), 
intent(inout) :: fs(n)
 
  192    real(kind=
rp), 
intent(in) :: b(n)
 
  193    real(kind=
rp), 
intent(in) :: dt, rho, bd(4)
 
  195    type(
field_t), 
pointer :: temp1, temp2
 
  196    integer :: temp_indices(2)
 
  201    do concurrent(i = 1:n)
 
  202       temp2%x(i,1,1,1) = s%x(i,1,1,1) * b(i) * bd(2)
 
  206       do concurrent(i = 1:n)
 
  207          temp1%x(i,1,1,1) = s_lag%lf(ilag-1)%x(i,1,1,1) * b(i) * bd(ilag+1)
 
  210       do concurrent(i = 1:n)
 
  211          temp2%x(i,1,1,1) = temp2%x(i,1,1,1) + temp1%x(i,1,1,1)
 
  215    do concurrent(i = 1:n)
 
  216       fs(i) = fs(i) + temp2%x(i,1,1,1) * (rho / dt)
 
 
  224    real(kind=
rp), 
intent(in) :: rho, dt
 
  225    integer, 
intent(in) :: n
 
  226    real(kind=
rp), 
intent(inout) :: bf_x(n), bf_y(n), bf_z(n)
 
  227    real(kind=
rp), 
intent(inout) :: phi_x(n), phi_y(n), phi_z(n)
 
  230    do concurrent(i = 1:n)
 
  231       bf_x(i) = bf_x(i) + phi_x(i) * (rho / dt)
 
  232       bf_y(i) = bf_y(i) + phi_y(i) * (rho / dt)
 
  233       bf_z(i) = bf_z(i) + phi_z(i) * (rho / dt)
 
 
subroutine rhs_maker_bdf_cpu(ulag, vlag, wlag, bfx, bfy, bfz, u, v, w, b, rho, dt, bd, nbd, n)