plan4.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
subroutine plan4()
  use kinds, only : dp
  use size_m, only : nid 
  use size_m, only : lx1, ly1, lz1, lelt
  use size_m, only : lx2, ly2, lz2, lelv
  use size_m, only : nx1, ny1, nz1, nelv
  use ctimer, only : icalld, tpres, npres, etime1, dnekclock
  use ctimer, only : np4misc, tp4misc
  use geom, only : binvm1, bm1, volvm1
  use helmholtz, only : hsolve, approx_space, init_approx_space
  use soln, only : vx, vy, vz, v1mask, v2mask, v3mask
  use soln, only : vtrans, pmask, vmult, pr
  use tstep, only : imesh, nmxh, tolhv
  use input, only : param
  implicit none

  integer, parameter :: ktot = lx1*ly1*lz1*lelt
  integer :: laxt, v_proj_size

  type(approx_space), save :: p_apx, vx_apx, vy_apx, vz_apx

  real(DP), allocatable :: RES1(:,:,:,:)
  real(DP), allocatable :: RES2(:,:,:,:)   
  real(DP), allocatable :: RES3(:,:,:,:)      
  real(DP), allocatable :: DV1 (:,:,:,:)   
  real(DP), allocatable :: DV2 (:,:,:,:)   
  real(DP), allocatable :: DV3 (:,:,:,:)   
  real(DP), allocatable :: RESPR (:,:,:,:)

  real(DP), allocatable :: h1(:,:,:,:), h2(:,:,:,:)
  real(DP), allocatable :: VEXT(:,:)
  REAL(DP), allocatable :: DPR(:,:,:,:)

  REAL(DP), allocatable :: DVC (:,:,:,:), DFC(:,:,:,:)
  REAL(DP) :: DIV1, DIV2, DIF1, DIF2, QTL1, QTL2
  real(DP) :: tolspl
  real(DP) :: etime
  real(DP), external :: glsum

  integer :: intype, ntot1 

  etime = dnekclock()
  np4misc = np4misc + 1
  v_proj_size = int(param(92))
  laxt = int(param(93))

  if (.not. allocated(p_apx%projectors)) then
    call init_approx_space(p_apx, laxt, ktot)
  endif
  if (.not. allocated(vx_apx%projectors)) then
    call init_approx_space(vx_apx, v_proj_size, ktot)
  endif
  if (.not. allocated(vy_apx%projectors)) then
    call init_approx_space(vy_apx, v_proj_size, ktot)
  endif
  if (.not. allocated(vz_apx%projectors)) then
    call init_approx_space(vz_apx, v_proj_size, ktot)
  endif

  intype = -1
  ntot1  = nx1*ny1*nz1*nelv

  ! add user defined divergence to qtl
!max  call add2 (qtl,usrdiv,ntot1)

  allocate(vext(lx1*ly1*lz1*lelv,3))

  ! Time-advance velocity with AB(k)
  CALL v_extrap(vext)

  ! compute explicit contributions (bf{x,y,z}) with BDF(k)
  etime = etime - dnekclock()
  CALL makef()
  etime = etime + dnekclock()

  ! store the velocity field for later
  CALL lagvel()

! split viscosity into explicit/implicit part
!    if (ifexplvis) call split_vis

! extrapolate velocity

  ! mask Dirichlet boundaries
  CALL bcdirvc(vx,vy,vz,v1mask,v2mask,v3mask)

!     first, compute pressure
#ifndef NOTIMER
  if (icalld == 0) tpres=0.0
  icalld=icalld+1
  npres=icalld
  etime1=dnekclock()
#endif

  allocate(respr(lx2,ly2,lz2,lelv))
  etime = etime - dnekclock()
  call crespsp(respr, vext)
  etime = etime + dnekclock()
  deallocate(vext)

  allocate(h1(lx1,ly1,lz1,lelv), h2(lx1,ly1,lz1,lelv))
  h1 = 1._dp / vtrans(:,:,:,:,1)
  h2 = 0._dp
  call ctolspl(tolspl,respr)

  allocate(dpr(lx2,ly2,lz2,lelv))
  etime = etime - dnekclock()
  call hsolve('PRES', dpr, respr, h1, h2, pmask, vmult &
  ,imesh,tolspl,nmxh,1 &
  ,p_apx,binvm1)
  etime = etime + dnekclock()
  deallocate(respr)
  pr = pr + dpr
  deallocate(dpr)
  call ortho(pr)
#ifndef NOTIMER
  tpres=tpres+(dnekclock()-etime1)
#endif

!     Compute velocity
  allocate(res1(lx1,ly1,lz1,lelv), &
           res2(lx1,ly1,lz1,lelv), &
           res3(lx1,ly1,lz1,lelv)  )

  etime = etime - dnekclock()
  call cresvsp(res1,res2,res3,h1,h2)
  etime = etime + dnekclock()

  allocate(dv1(lx1,ly1,lz1,lelv), &
           dv2(lx1,ly1,lz1,lelv), &
           dv3(lx1,ly1,lz1,lelv)  )

  etime = etime - dnekclock()
  call hsolve('VELX', dv1, res1, h1, h2, v1mask, vmult, imesh, tolhv, nmxh, 1, &
              vx_apx, binvm1)
  call hsolve('VELY', dv2, res2, h1, h2, v2mask, vmult, imesh, tolhv, nmxh, 2, &
              vy_apx, binvm1)
  call hsolve('VELZ', dv3, res3, h1, h2, v3mask, vmult, imesh, tolhv, nmxh, 3, &
              vz_apx, binvm1)
  etime = etime + dnekclock()
  deallocate(res1, res2, res3, h1, h2)

  vx = vx + dv1; vy = vy + dv2; vz = vz + dv3

!    if (ifexplvis) call redo_split_vis

! Below is just for diagnostics...

!     Calculate Divergence norms of new VX,VY,VZ
  allocate(dvc(lx1,ly1,lz1,lelv), dfc(lx1,ly1,lz1,lelv))
  CALL opdiv(dvc,vx,vy,vz)
  CALL dssum(dvc)
  dvc = dvc * binvm1

  dv1 = dvc * bm1
  div1 = glsum(dv1,ntot1)/volvm1

  dv2 = dvc * dvc * bm1
  div2 = glsum(dv2,ntot1)/volvm1
  div2 = sqrt(div2)
!     Calculate Divergence difference norms
  dfc = dvc! - qtl
  dv1 = dfc * bm1
  dif1 = glsum(dv1,ntot1)/volvm1
    
  dv2 = dfc * dfc * bm1
  dif2 = glsum(dv2,ntot1)/volvm1
  dif2 = sqrt(dif2)

  dv1 = 0._dp !qtl * bm1
  qtl1 = glsum(dv1,ntot1)/volvm1
    
  dv2 = 0._dp !qtl * qtl * bm1
  qtl2 = glsum(dv2,ntot1)/volvm1
  qtl2 = sqrt(qtl2)

  IF (nid == 0) THEN
      WRITE(6,'(15X,A,1p2e13.4)') &
      'L1/L2 DIV(V)    :',div1,div2
      WRITE(6,'(15X,A,1p2e13.4)') &
      'L1/L2 QTL       :',qtl1,qtl2
      WRITE(6,'(15X,A,1p2e13.4)') &
      'L1/L2 DIV(V)-QTL:',dif1,dif2
      IF (dif2 > 0.1) WRITE(6,'(15X,A)') &
      'WARNING: DIV(V)-QTL too large!'
  ENDIF 
  tp4misc = tp4misc + (dnekclock() - etime)
   
  return
end subroutine plan4

!-----------------------------------------------------------------------
subroutine crespsp (respr, vext)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lx2, ly2, lz2, lelv
  use size_m, only : nx1, ny1, nz1, nelv, ndim, nx2, ny2, nz2
  use geom, only : rxm2, sxm2, txm2, rym2, sym2, tym2, rzm2, szm2, tzm2
  use geom, only : unx, uny, unz, area
  use input, only : ifaxis, if3d, cbc
  use geom, only : bm1, binvm1, jacmi
  use soln, only : bfx, bfy, bfz, pr
  use soln, only : vx, vy, vz
  use soln, only : vtrans, vdiff, vmult!, qtl
  use soln, only : v1mask, v2mask, v3mask
  use tstep, only : imesh, bd, dt, ifield
  use mesh, only : if_ortho
  use dxyz, only : dxtm12, dym12, dzm12
  use ctimer, only : ncrespsp, tcrespsp, dnekclock
  use parallel, only : nid
  implicit none

  REAL(DP), intent(out) :: RESPR (lx2,ly2,lz2,lelv)
  real(DP), intent(in)  :: VEXT  (lx1*ly1*lz1*lelv,3)

  real(DP), allocatable, dimension(:,:,:,:) :: TA1, TA2, TA3
  real(DP), allocatable, dimension(:,:,:,:) :: WA1, WA2, WA3
  real(DP), allocatable, dimension(:,:,:,:) :: W1,  W2

  CHARACTER(3) :: CB
 
  integer :: nxyz1, ntot1, nfaces       
  integer :: n, ifc, iel, e, iz
  integer :: nyz2, nxy1
  real(DP) :: scale, dtbd
  real(DP) :: etime
  real(DP), allocatable :: tmp1(:,:,:), tmp2(:,:,:), tmp3(:,:,:)
  logical :: any_face

  ncrespsp = ncrespsp + 1
  etime = dnekclock()

  allocate(ta1(lx1,ly1,lz1,lelv) &
  , ta2(lx1,ly1,lz1,lelv) &
  , ta3(lx1,ly1,lz1,lelv) )
  allocate(w1(lx1,ly1,lz1,lelv) &
  , w2(lx1,ly1,lz1,lelv) )

 
  nxyz1  = nx1*ny1*nz1
  ntot1  = nxyz1*nelv
  nfaces = 2*ndim

!   -mu*curl(curl(v))
  call op_curl(ta1,ta2,ta3,vext(1,1),vext(1,2),vext(1,3), &
   .true. ,w1,w2)
  if(ifaxis) then
!max      CALL COL2 (TA2, OMASK,NTOT1)
!max      CALL COL2 (TA3, OMASK,NTOT1)
  endif

  allocate( wa1(lx1,ly1,lz1,lelv) &
  , wa2(lx1,ly1,lz1,lelv) &
  , wa3(lx1,ly1,lz1,lelv) )
  call op_curl(wa1,wa2,wa3,ta1,ta2,ta3, .true. ,w1,w2)
  deallocate(w1, w2)

  allocate(tmp1(lx1, ly1, lz1), tmp2(lx1,ly1,lz1), tmp3(lx1,ly1,lz1))

  if(ifaxis) then
!max      CALL COL2  (WA2, OMASK,NTOT1)
!max      CALL COL2  (WA3, OMASK,NTOT1)
  endif
  !wa1 = wa1 * bm1; wa2 = wa2 * bm1; wa3 = wa3 * bm1

  !call opgrad   (ta1,ta2,ta3,QTL)
  !ta1 = 0._dp; ta2 = 0._dp; ta3 = 0._dp
  if(ifaxis) then
!max      CALL COL2  (ta2, OMASK,ntot1)
!max      CALL COL2  (ta3, OMASK,ntot1)
  endif

!   add old pressure term because we solve for delta p
  ta2 = 0._dp
  ta1 = 1._dp / vtrans(:,:,:,:,1)

  etime = etime - dnekclock()
  CALL axhelm(respr,pr,ta1,ta2,imesh,1)
  etime = etime + dnekclock()

  !   add explicit (NONLINEAR) terms
  n = nx1*ny1*nz1*nelv
  do iel = 1, nelv
    !tmp1 = vdiff(:,:,:,iel,1) *ta1(:,:,:,iel)
    tmp1 = bm1(:,:,:,iel) * vdiff(:,:,:,iel,1) *ta1(:,:,:,iel)
    ta3(:,:,:,iel) = binvm1(:,:,:,iel) * (bfz(:,:,:,iel)*ta1(:,:,:,iel)-tmp1*wa3(:,:,:,iel))
    ta2(:,:,:,iel) = binvm1(:,:,:,iel) * (bfy(:,:,:,iel)*ta1(:,:,:,iel)-tmp1*wa2(:,:,:,iel))
    ta1(:,:,:,iel) = binvm1(:,:,:,iel) * (bfx(:,:,:,iel)*ta1(:,:,:,iel)-tmp1*wa1(:,:,:,iel))
  enddo

  call opdssum(ta1,ta2,ta3)

  if (if3d) then

    if (if_ortho) then
      !deallocate(wa1,wa2,wa3)
      nyz2  = ny2*nz2
      nxy1  = nx1*ny1

     do e = 1, nelv
        tmp1 = jacmi(:,:,:,e) * bm1(:,:,:,e)
        ! X 
        tmp2 = ta1(:,:,:,e) * rxm2(:,:,:,e) * tmp1
        call mxm(dxtm12,nx1,tmp2,nx2,tmp3,nyz2)
        respr(:,:,:,e) = - respr(:,:,:,e) + tmp3
        ! Y 
        tmp2 = ta2(:,:,:,e) * sym2(:,:,:,e) * tmp1 
        do iz=1,nz2
            call mxm(tmp2(:,:,iz),nx1,dym12,ny2,tmp3(:,:,iz),ny1)
        enddo
        respr(:,:,:,e) =   respr(:,:,:,e) + tmp3
        ! Z
        tmp2 = ta3(:,:,:,e) * tzm2(:,:,:,e) * tmp1 
        call mxm(tmp2,nxy1,dzm12,nz2,tmp3,nz1)
        respr(:,:,:,e) =   respr(:,:,:,e) + tmp3
      enddo

    else
      call cdtp(wa1,ta1,rxm2,sxm2,txm2,1)
      call cdtp(wa2,ta2,rym2,sym2,tym2,1)
      call cdtp(wa3,ta3,rzm2,szm2,tzm2,1)
      respr = -respr + wa1 + wa2 + wa3
      !deallocate(wa1,wa2,wa3)
    endif
  else
      call cdtp(wa1,ta1,rxm2,sxm2,txm2,1)
      call cdtp(wa2,ta2,rym2,sym2,tym2,1)
      respr = -respr + wa1 + wa2 
  endif
  deallocate(wa1,wa2,wa3)

!   add thermal divergence
  dtbd = bd(1)/dt
  respr = respr !+ qtl * bm1 * dtbd
   
!   surface terms
  DO iel=1,nelv
      DO ifc=1,nfaces
          cb = cbc(ifc,iel,ifield)
          any_face = .false.
          IF (cb(1:1) == 'V' .OR. cb(1:1) == 'v') THEN
              any_face = .true.
              tmp1 = 0._dp; tmp2 = 0._dp
              IF (ndim == 3) tmp3 = 0._dp
              CALL faccl3 &
              (tmp1,vx(1,1,1,iel),unx(1,1,ifc,iel),ifc)
              CALL faccl3 &
              (tmp2,vy(1,1,1,iel),uny(1,1,ifc,iel),ifc)
              IF (ndim == 3) &
              CALL faccl3 &
              (tmp3,vz(1,1,1,iel),unz(1,1,ifc,iel),ifc)
          else if (cb(1:3) == 'SYM') then
            any_face = .true.
            tmp1 = 0._dp; tmp2 = 0._dp
            IF (ndim == 3) tmp3 = 0._dp

            CALL faccl3(tmp1,ta1(:,:,:,iel),unx(1,1,ifc,iel),ifc)
            CALL faccl3(tmp2,ta2(:,:,:,iel),uny(1,1,ifc,iel),ifc)
            IF (ndim == 3) then
              CALL faccl3(tmp3,ta3(:,:,:,iel),unz(1,1,ifc,iel),ifc)
            endif
          endif

          if (any_face) then
            IF (ndim == 3) then 
              tmp1 = tmp1 + tmp2 + tmp3
            else
              tmp1 = tmp1 + tmp2
            endif
            CALL faccl2(tmp1,area(1,1,ifc,iel),ifc)
          endif

          if (cb(1:1) == 'V' .OR. cb(1:1) == 'v') then
            respr(:,:,:,iel) = respr(:,:,:,iel) - dtbd*tmp1
          else if (cb(1:3) == 'SYM') then
            respr(:,:,:,iel) = respr(:,:,:,iel) - tmp1 
          endif
      END DO
      !call dssum(ta1) ! maybe this should be here for SYM? (maxhutch)
  END DO
  deallocate(ta1, ta2, ta3)

!   Assure that the residual is orthogonal to (1,1,...,1)T
!   (only if all Dirichlet b.c.)
  CALL ortho(respr)
  tcrespsp = tcrespsp + (dnekclock() - etime)

  return
end subroutine crespsp

!----------------------------------------------------------------------
subroutine cresvsp (resv1,resv2,resv3,h1,h2)
  use kinds, only : dp
  use size_m, only : nx1, ny1, nz1, nelv
  use size_m, only : lx1, ly1, lz1, lelv
  use input, only : ifaxis
  use soln, only : vx, vy, vz, pr, bfx, bfy, bfz !, qtl
  use ctimer, only : ncresvsp, tcresvsp, dnekclock
  implicit none

  real(DP), intent(out) :: resv1(lx1,ly1,lz1,lelv) 
  real(DP), intent(out) :: resv2(lx1,ly1,lz1,lelv) 
  real(DP), intent(out) :: resv3(lx1,ly1,lz1,lelv) 
  real(DP), intent(out)  :: h1   (lx1,ly1,lz1,lelv) 
  real(DP), intent(out)  :: h2   (lx1,ly1,lz1,lelv)

  real(DP), allocatable, dimension(:,:,:,:) :: TA1, TA2, TA3, TA4

  integer :: ntot, intype
  real(DP) :: scale
  real(DP) :: etime

  etime = dnekclock()
  ncresvsp = ncresvsp + 1

  ntot = nx1*ny1*nz1*nelv
  intype = -1

  CALL sethlm(h1,h2,intype)

  etime = etime - dnekclock()
  ! just three axhelm calls
  CALL ophx(resv1,resv2,resv3,vx,vy,vz,h1,h2)
  etime = etime + dnekclock()

  scale = -1./3.
  allocate(ta1(lx1,ly1,lz1,lelv) &
  ,             ta2(lx1,ly1,lz1,lelv) &
  ,             ta3(lx1,ly1,lz1,lelv) &
  ,             ta4(lx1,ly1,lz1,lelv) )
  ta4 = pr !+ scale*(vdiff(:,:,:,:,1) * qtl)

  call opgrad(ta1,ta2,ta3,ta4)
  deallocate(ta4)

  if(ifaxis) then
!max      CALL COL2 (TA2, OMASK,NTOT)
!max      CALL COL2 (TA3, OMASK,NTOT)
  endif

  resv1 = -resv1 + bfx - ta1
  resv2 = -resv2 + bfy - ta2
  resv3 = -resv3 + bfz - ta3

  tcresvsp = tcresvsp + (dnekclock() - etime)

  return
end subroutine cresvsp

!-----------------------------------------------------------------------
subroutine op_curl(w1,w2,w3,u1,u2,u3,ifavg,work1,work2)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lelv, nx1, ny1, nz1, nelv
  use geom, only : rxm1, rym1, rzm1, sxm1, sym1, szm1, txm1, tym1, tzm1
  use dxyz, only : dztm1, dytm1, dxm1
  use geom, only : jacm1, bm1, binvm1, jacmi
  use input, only : if3d, ifaxis, ifcyclic
  use tstep, only : ifield
  use mesh, only : if_ortho
  implicit none

  real(DP), intent(out) :: w1(lx1,ly1,lz1,lelv) !>
  real(DP), intent(out) :: w2(lx1,ly1,lz1,lelv) !>
  real(DP), intent(out) :: w3(lx1,ly1,lz1,lelv) !>
  real(DP), intent(in)  :: u1(lx1,ly1,lz1,lelv) !>
  real(DP), intent(in)  :: u2(lx1,ly1,lz1,lelv) !>
  real(DP), intent(in)  :: u3(lx1,ly1,lz1,lelv) !>
  real(DP), intent(out) :: work1(lx1,ly1,lz1,lelv) !>
  real(DP), intent(out) :: work2(lx1,ly1,lz1,lelv) !>
  logical, intent(in)   :: ifavg !>

  integer :: ntot, nxyz, ifielt, nxy1, nyz1, iel, iz
  real(DP), allocatable :: tmp1(:,:,:), tmp2(:,:,:)

  allocate(tmp1(nx1,ny1,nz1), tmp2(nx1,ny1,nz1))

  ntot  = nx1*ny1*nz1*nelv
  nxyz  = nx1*ny1*nz1
  nxy1  = nx1*ny1
  nyz1  = ny1*nz1

!   work1=dw/dy ; work2=dv/dz
  if (if3d) then
    if (if_ortho) then
      do iel = 1, nelv
        do iz = 1, nz1
          CALL mxm(u3(1,1,iz,iel),nx1,dytm1,ny1,tmp1(1,1,iz),ny1)
        enddo
        CALL mxm(u2(1,1,1,iel),nxy1,dztm1,nz1,tmp2,nz1) 
        w1(:,:,:,iel) = (tmp1*sym1(:,:,:,iel) - tmp2*tzm1(:,:,:,iel)) * jacmi(:,:,:,iel)
      enddo
    else
      call dudxyz(work1,u3,rym1,sym1,tym1,jacm1,1,2)
      call dudxyz(work2,u2,rzm1,szm1,tzm1,jacm1,1,3)
      w1 = work1(:,:,:,1:nelv) - work2(:,:,:,1:nelv) 
    endif
  else
      call dudxyz(work1,u3,rym1,sym1,tym1,jacm1,1,2)
      call copy(w1,work1,ntot)

      if(ifaxis) then
        write(*,*) "Oops: ifaxis"
#if 0
          call copy(ta,u3,ntot)
          do iel = 1,nelv
              if(ifrzer(iel)) then
                  call rzero(ta(1,1,1,iel),nx1)
                  call mxm(ta(1,1,1,iel),nx1,datm1,ny1,duax,1)
                  call copy(ta(1,1,1,iel),duax,nx1)
              endif
              call col2(ta(1,1,1,iel),yinvm1(1,1,1,iel),nxyz)
          enddo
          call add2(w1,ta,ntot)
#endif
      endif
  endif

!   work1=du/dz ; work2=dw/dx
  if (if3d) then
      if (if_ortho) then
        do iel = 1, nelv
          CALL mxm(u1(1,1,1,iel),nxy1,dztm1,nz1,tmp1,nz1) 
          CALL mxm(dxm1,nx1,u3(1,1,1,iel),nx1,tmp2,nyz1)
          w2(:,:,:,iel) = (tmp1*tzm1(:,:,:,iel) - tmp2*rxm1(:,:,:,iel)) * jacmi(:,:,:,iel)
        enddo
      else
        call dudxyz(work1,u1,rzm1,szm1,tzm1,jacm1,1,3)
        call dudxyz(work2,u3,rxm1,sxm1,txm1,jacm1,1,1)
        w2 = work1(:,:,:,1:nelv) - work2(:,:,:,1:nelv)
      endif
  else
      work1(:,:,:,1:nelv) = 0._dp
      call dudxyz(work2,u3,rxm1,sxm1,txm1,jacm1,1,1)
      w2 = work1(:,:,:,1:nelv) - work2(:,:,:,1:nelv)
  endif

!   work1=dv/dx ; work2=du/dy
  if (if_ortho) then
    do iel = 1, nelv
      CALL mxm(dxm1,nx1,u2(1,1,1,iel),nx1,tmp1,nyz1)
      do iz = 1, nz1
        CALL mxm(u1(1,1,iz,iel),nx1,dytm1,ny1,tmp2(1,1,iz),ny1)
      enddo
      w3(:,:,:,iel) = (tmp1*rxm1(:,:,:,iel) - tmp2*sym1(:,:,:,iel)) * jacmi(:,:,:,iel)
    enddo
  else
    call dudxyz(work1,u2,rxm1,sxm1,txm1,jacm1,1,1)
    call dudxyz(work2,u1,rym1,sym1,tym1,jacm1,1,2)
    w3 = work1(:,:,:,1:nelv) - work2(:,:,:,1:nelv)
  endif
!  Avg at bndry

!   if (ifavg) then
  if (ifavg .AND. .NOT. ifcyclic) then

      ifielt = ifield
      ifield = 1
             
      w1 = w1 * bm1; w2 = w2 * bm1; w3 = w3 * bm1
      call opdssum(w1,w2,w3)
      w1 = w1 * binvm1; w2 = w2 * binvm1; w3 = w3 * binvm1

      ifield = ifielt

  endif

  return
end subroutine op_curl

!-----------------------------------------------------------------------
subroutine v_extrap(vext)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lelv
  use input, only : if3d
  use soln, only : vx, vy, vz, vxlag, vylag, vzlag
  use tstep, only : ab, nab
  implicit none

  real(DP), intent(out) :: vext(lx1,ly1,lz1,lelv,3) 
  real(DP) :: AB0, AB1, AB2

  ab0 = ab(1)
  ab1 = ab(2)
  ab2 = ab(3)

  if(nab == 3) then
             vext(:,:,:,:,1) = ab0*vx + ab1*vxlag(:,:,:,:,1) + ab2*vxlag(:,:,:,:,2)
             vext(:,:,:,:,2) = ab0*vy + ab1*vylag(:,:,:,:,1) + ab2*vylag(:,:,:,:,2)
    if(if3d) vext(:,:,:,:,3) = ab0*vz + ab1*vzlag(:,:,:,:,1) + ab2*vzlag(:,:,:,:,2)
  else
             vext(:,:,:,:,1) = ab0*vx + ab1*vxlag(:,:,:,:,1)
             vext(:,:,:,:,2) = ab0*vy + ab1*vylag(:,:,:,:,1)
    if(if3d) vext(:,:,:,:,3) = ab0*vz + ab1*vzlag(:,:,:,:,1)
  endif

  return
end subroutine v_extrap

!-----------------------------------------------------------------------