navier1.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
subroutine ctolspl (tolspl,respr)
  use kinds, only : dp
  use size_m, only : lelv, lx2, ly2, lz2
  use size_m, only : nx1, ny1, nz1, nelv, nid
  use geom, only : binvm1, volvm1
  use tstep, only : dt, tolpdf, tolps, prelax
  use soln, only : pmask, vmult
  implicit none

  real(DP) :: tolspl
  REAL(DP) :: RESPR (lx2,ly2,lz2,lelv)

  real(DP), allocatable :: WORK(:,:,:,:)
  integer :: ntot1
  real(DP) :: rinit, tolmin, tolold
  real(DP), external :: glsc23

  allocate(work(lx2,ly2,lz2,lelv))
  ntot1 = nx1*ny1*nz1*nelv
  work = respr
  call dssum(work)
  work = work * pmask
  rinit = glsc23(work, binvm1, vmult, ntot1)
  rinit  = sqrt(rinit/volvm1)
  IF (tolpdf > 0.) THEN
      tolspl = tolpdf
      tolmin = tolpdf
  ELSE
      tolspl = tolps/dt
      tolmin = rinit*prelax
  ENDIF
  IF (tolspl < tolmin) THEN
      tolold = tolspl
      tolspl = tolmin
      IF (nid == 0) &
      WRITE (6,*) 'Relax the pressure tolerance ',tolspl,tolold
  ENDIF

  return
end subroutine ctolspl

!------------------------------------------------------------------------
subroutine ortho (respr)
  use kinds, only : dp, i8, i4
  use size_m, only : lx2, ly2, lz2, lelv, nx2, ny2, nz2, nelv
  use geom, only : ifvcor, ifbcor
  use input, only : ifldmhd
  use parallel, only : nelgv
  use tstep, only : ifield
  implicit none

  real(DP) :: respr (lx2,ly2,lz2,lelv)
  integer(i8) :: ntotg,nxyz2
  integer :: ntot
  real(DP) :: rlam
  real(DP), external :: glsum

  nxyz2 = nx2*ny2*nz2
  ntot  = int(nxyz2*nelv, kind=i4)
  ntotg = nxyz2*nelgv

  if (ifield == 1) then
      if (ifvcor) then
          rlam  = glsum(respr,ntot)/ntotg
          respr = respr - rlam
      endif
  elseif (ifield == ifldmhd) then
      if (ifbcor) then
          rlam = glsum(respr,ntot)/ntotg
          respr = respr - rlam
      endif
  else
      call exitti('ortho: unaccounted ifield = $',ifield)
  endif

  return
end subroutine ortho

!---------------------------------------------------------------------
!---------------------------------------------------------------------
subroutine opgrad (out1,out2,out3,inp)
  use kinds, only : dp
  use size_m, only : nx2, ny2, nz2, ndim, nelv
  use size_m, only : lx1, ly1, lz1, lx2, ly2, lz2, nx1, ny1, nz1
  use geom, only : rxm2, rym2, rzm2, sxm2, sym2, szm2, txm2, tym2, tzm2
  use geom, only : bm1, jacmi
  use input, only : ifsplit, ifaxis
  use mesh, only : if_ortho
  use dxyz, only : dxm12, dytm12, dztm12
  use ixyz, only : ixm12, iytm12, iztm12
  implicit none

  REAL(DP) :: OUT1 (lx2,ly2,lz2,nelv)
  REAL(DP) :: OUT2 (lx2,ly2,lz2,nelv)
  REAL(DP) :: OUT3 (lx2,ly2,lz2,nelv)
  REAL(DP) :: INP  (lx1,ly1,lz1,nelv)

  real(DP) ::  ta1 (lx1*ly1*lz1) &
  ,             ta2 (lx1*ly1*lz1) &
  ,             ta3 (lx1*ly1*lz1)

  integer :: iflg, ntot2, i1, i2, e, iz, n1, n2, nxy2, nyz1
  iflg = 0

  if (ifsplit .AND. .NOT. ifaxis) then
      call wgradm1(out1,out2,out3,inp,nelv) ! weak grad on FLUID mesh
      return
  endif

  ntot2 = nx2*ny2*nz2*nelv
  if (if_ortho) then
    nxy2  = nx2*ny2
    n1    = nx2*ny1
    n2    = nx2*ny2
    nyz1  = ny1*nz1

    do e=1,nelv
      call mxm(dxm12,nx2,inp(:,:,:,e),nx1,ta1,nyz1)
      i1=1
      i2=1
      do iz=1,nz1
          call mxm(ta1(i1),nx2,iytm12,ny1,ta2(i2),ny2)
          i1=i1+n1
          i2=i2+n2
      enddo
      call mxm(ta2,nxy2,iztm12,nz1,out1(:,:,:,e),nz2)    
 
      call mxm(ixm12,nx2,inp(:,:,:,e),nx1,ta3,nyz1) ! reuse ta3 below
      i1=1
      i2=1
      do iz=1,nz1
          call mxm(ta3(i1),nx2,dytm12,ny1,ta2(i2),ny2)
          i1=i1+n1
          i2=i2+n2
      enddo
      call mxm(ta2,nxy2,iztm12,nz1,out2(:,:,:,e),nz2)

      ! re-use ta3 from above
      i1=1
      i2=1
      do iz=1,nz1
          call mxm(ta3(i1),nx2,iytm12,ny1,ta2(i2),ny2)
          i1=i1+n1
          i2=i2+n2
      enddo
      call mxm(ta2,nxy2,dztm12,nz1,out3(:,:,:,e),nz2)
    enddo

    out1 = out1 * rxm2 * bm1 * jacmi
    out2 = out2 * sym2 * bm1 * jacmi 
    out3 = out3 * tzm2 * bm1 * jacmi

  else
    CALL multd(out1,inp,rxm2,sxm2,txm2,1,iflg)
    CALL multd(out2,inp,rym2,sym2,tym2,2,iflg)
    IF (ndim == 3) &
    CALL multd(out3,inp,rzm2,szm2,tzm2,3,iflg)
  endif

  return
end subroutine opgrad

!-------------------------------------------------------------
!-------------------------------------------------------------
subroutine cdtp (dtx,x,rm2,sm2,tm2,isd)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lx2, ly2, lz2, lelv
  use size_m, only : nx1, ny1, nz1, nx2, ny2, nz2, nelv, ndim
  use ctimer, only : icalld, tcdtp, ncdtp, etime1, dnekclock
  use dxyz, only : dym12, dam12, dcm12, dxtm12, dzm12
  use geom, only : ifrzer, jacm2, ym2, jacm1
  use input, only : ifaxis, ifsplit
  use ixyz, only : iym12, iam12, icm12
  use geom, only : bm1, bm2
  use wz_m, only : w3m2, w2am2, w2cm2
  implicit none

  integer :: isd
  real(DP) :: dtx  (lx1,ly1,lz1,lelv)
  real(DP) :: x    (lx2,ly2,lz2,lelv)
  real(DP) :: rm2  (lx2,ly2,lz2,lelv)
  real(DP) :: sm2  (lx2,ly2,lz2,lelv)
  real(DP) :: tm2  (lx2,ly2,lz2,lelv)

  real(DP) ::  wx  (lx1,ly1,lz1) &
  ,             ta1 (lx1,ly1,lz1) &
  ,             ta2 (lx1,ly1,lz1)

  integer :: e
  integer :: nxyz1, nxyz2, nxy1, nyz2, n1, n2, ny12, i1, i2, iz

#ifndef NOTIMER
  if (icalld == 0) tcdtp=0.0
  icalld=icalld+1
  ncdtp=icalld
  etime1=dnekclock()
#endif

  nxyz1 = nx1*ny1*nz1
  nxyz2 = nx2*ny2*nz2
  nyz2  = ny2*nz2
  nxy1  = nx1*ny1

  n1    = nx1*ny1
  n2    = nx1*ny2

  do e=1,nelv

  !       Use the appropriate derivative- and interpolation operator in
  !       the y-direction (= radial direction if axisymmetric).
      if (ifaxis) then
          ny12   = ny1*ny2
          if (ifrzer(e)) then
              call copy(iym12,iam12,ny12)
              call copy(dym12,dam12,ny12)
              call copy(w3m2,w2am2,nxyz2)
          else
              call copy(iym12,icm12,ny12)
              call copy(dym12,dcm12,ny12)
              call copy(w3m2,w2cm2,nxyz2)
          endif
      endif
  
  !      Collocate with weights
  
      if(ifsplit) then
          wx = bm1(:,:,:,e) * x(:,:,:,e) / jacm1(:,:,:,e)
      else
          if (ifaxis) then
              if (ifrzer(e)) then
                  wx = x(:,:,:,e) * bm2(:,:,:,e) / jacm2(:,:,:,e)
              else
                  wx = w3m2 * x(:,:,:,e) * ym2(:,:,:,e)
              endif
          else
            wx = w3m2 * x(:,:,:,e) 
          endif
      endif
  
      if (ndim == 2) then
          write(*,*) "Whoops! cdtp"
#if 0
          if ( .NOT. ifdfrm(e) .AND. ifalgn(e)) then
          
              if (      ifrsxy(e) .AND. isd == 1  .OR. &
               .NOT. ifrsxy(e) .AND. isd == 2) then
              
                  call col3(ta1,wx,rm2(1,e),nxyz2)
                  call mxm(dxtm12,nx1,ta1,nx2,ta2,nyz2)
                  call mxm(ta2,nx1,iym12,ny2,dtx(1,e),ny1)
              else
                  call col3(ta1,wx,sm2(1,e),nxyz2)
                  call mxm(ixtm12,nx1,ta1,nx2,ta2,nyz2)
                  call mxm(ta2,nx1,dym12,ny2,dtx(1,e),ny1)
              endif
          else
              call col3(ta1,wx,rm2(1,e),nxyz2)
              call mxm(dxtm12,nx1,ta1,nx2,ta2,nyz2)
              call mxm(ta2,nx1,iym12,ny2,dtx(1,e),ny1)

              call col3(ta1,wx,sm2(1,e),nxyz2)
              call mxm(ixtm12,nx1,ta1,nx2,ta2,nyz2)
              call mxm(ta2,nx1,dym12,ny2,ta1,ny1)

              call add2(dtx(1,e),ta1,nxyz1)
          endif
#endif
      else
          if (ifsplit) then
              ta1 = wx * rm2(:,:,:,e)
              call mxm(dxtm12,nx1,ta1,nx2,dtx(:,:,:,e),nyz2)
              ta1 = wx * sm2(:,:,:,e) 
              i1 = 1
              i2 = 1
              do iz=1,nz2
                  call mxm(ta1(:,:,iz),nx1,dym12,ny2,ta2(:,:,iz),ny1)
                  i1 = i1 + n1
                  i2 = i2 + n2
              enddo
              dtx(:,:,:,e) = dtx(:,:,:,e) + ta2
              ta1 = wx * tm2(:,:,:,e)
              call mxm(ta1,nxy1,dzm12,nz2,ta2,nz1)
              dtx(:,:,:,e) = dtx(:,:,:,e) + ta2
          else
              write(*,*) "Whoops! cdtp"
#if 0
              call col3(ta1,wx,rm2(1,e),nxyz2)
              call mxm(dxtm12,nx1,ta1,nx2,ta2,nyz2)
              i1 = 1
              i2 = 1
              do iz=1,nz2
                  call mxm(ta2(i2),nx1,iym12,ny2,ta1(i1),ny1)
                  i1 = i1 + n1
                  i2 = i2 + n2
              enddo
              call mxm(ta1,nxy1,izm12,nz2,dtx(1,e),nz1)

              call col3(ta1,wx,sm2(1,e),nxyz2)
              call mxm(ixtm12,nx1,ta1,nx2,ta2,nyz2)
              i1 = 1
              i2 = 1
              do iz=1,nz2
                  call mxm(ta2(i2),nx1,dym12,ny2,ta1(i1),ny1)
                  i1 = i1 + n1
                  i2 = i2 + n2
              enddo
              call mxm(ta1,nxy1,izm12,nz2,ta2,nz1)
              call add2(dtx(1,e),ta2,nxyz1)

              call col3(ta1,wx,tm2(1,e),nxyz2)
              call mxm(ixtm12,nx1,ta1,nx2,ta2,nyz2)
              i1 = 1
              i2 = 1
              do iz=1,nz2
                  call mxm(ta2(i2),nx1,iym12,ny2,ta1(i1),ny1)
                  i1 = i1 + n1
                  i2 = i2 + n2
              enddo
              call mxm(ta1,nxy1,dzm12,nz2,ta2,nz1)
              call add2(dtx(1,e),ta2,nxyz1)
#endif
          endif

      endif
  
  !     If axisymmetric, add an extra diagonal term in the radial
  !     direction (only if solving the momentum equations and ISD=2)
  !     NOTE: NZ1=NZ2=1
  if (ifaxis) write(*,*) "Oops: ifaxis" 
#if 0 
      if(ifsplit) then

          if (ifaxis .AND. (isd == 4)) then
              call copy(ta1,x(1,e),nxyz1)
              if (ifrzer(e)) THEN
                  call rzero(ta1, nx1)
                  call mxm(x(1,e),nx1,datm1,ny1,duax,1)
                  call copy(ta1,duax,nx1)
              endif
              call col2(ta1,baxm1(1,1,1,e),nxyz1)
              call add2(dtx(1,e),ta1,nxyz1)
          endif

      else

          if (ifaxis .AND. (isd == 2)) then
              call col3(ta1,x(1,e),bm2(1,1,1,e),nxyz2)
              ta = ta / ym2(:,:,:,e)
              call mxm(ixtm12,nx1,ta1,nx2,ta2,ny2)
              call mxm(ta2,nx1,iym12,ny2,ta1,ny1)
              call add2(dtx(1,e),ta1,nxyz1)
          endif

      endif
#endif
  enddo

#ifndef NOTIMER
  tcdtp=tcdtp+(dnekclock()-etime1)
#endif
  return
end subroutine cdtp

!---------------------------------------------------------------------
!---------------------------------------------------------------------
subroutine multd (dx,x,rm2,sm2,tm2,isd,iflg)
  use kinds, only : dp
  use size_m, only : lx2, ly2, lz2, lx1, ly1, lz1, lelv
  use size_m, only : nx1, ny1, nz1, nx2, ny2, nz2, nelv, ndim
  use ctimer, only : icalld, tmltd, nmltd, etime1, dnekclock
  use dxyz, only : dytm12, datm12, dctm12, dxm12, dztm12
  use geom, only : ifrzer, jacm1
  use input, only : ifaxis, ifsplit
  use ixyz, only : iytm12, iatm12, ictm12, iztm12, ixm12
  use geom, only : bm1
  use wz_m, only : w3m2, w2am2, w2cm2
  implicit none

  integer :: isd, iflg
  real(DP) ::           dx   (lx2,ly2,lz2,lelv)
  real(DP) ::           x    (lx1,ly1,lz1,lelv)
  real(DP) ::           rm2  (lx2,ly2,lz2,lelv)
  real(DP) ::           sm2  (lx2,ly2,lz2,lelv)
  real(DP) ::           tm2  (lx2,ly2,lz2,lelv)

  real(DP) ::  ta1 (lx1*ly1*lz1) &
  ,             ta2 (lx1*ly1*lz1) &
  ,             ta3 (lx1*ly1*lz1)

  integer :: e
  integer :: nxyz1, nxy2, nxyz2, n1, n2, ny12, i1, i2, iz, nyz1

#ifndef NOTIMER
  if (icalld == 0) tmltd=0.0
  icalld=icalld+1
  nmltd=icalld
  etime1=dnekclock()
#endif

  nyz1  = ny1*nz1
  nxy2  = nx2*ny2
  nxyz1 = nx1*ny1*nz1
  nxyz2 = nx2*ny2*nz2

  n1    = nx2*ny1
  n2    = nx2*ny2

  do e=1,nelv

  !        Use the appropriate derivative- and interpolation operator in
  !        the y-direction (= radial direction if axisymmetric).
      if (ifaxis) then
          ny12   = ny1*ny2
          if (ifrzer(e)) then
              call copy(iytm12,iatm12,ny12)
              call copy(dytm12,datm12,ny12)
              call copy(w3m2,w2am2,nxyz2)
          else
              call copy(iytm12,ictm12,ny12)
              call copy(dytm12,dctm12,ny12)
              call copy(w3m2,w2cm2,nxyz2)
          endif
      endif

      if (ndim == 2) then
          write(*,*) "Whoops! multd"
#if 0
          if ( .NOT. ifdfrm(e) .AND. ifalgn(e)) then
          
              if (      ifrsxy(e) .AND. isd == 1  .OR. &
               .NOT. ifrsxy(e) .AND. isd == 2) then
                  call mxm(dxm12,nx2,x(1,e),nx1,ta1,nyz1)
                  call mxm(ta1,nx2,iytm12,ny1,dx(1,e),ny2)
                  call col2(dx(1,e),rm2(1,e),nxyz2)
              else
                  call mxm(ixm12,nx2,x(1,e),nx1,ta1,nyz1)
                  call mxm(ta1,nx2,dytm12,ny1,dx(1,e),ny2)
                  call col2(dx(1,e),sm2(1,e),nxyz2)
              endif
          else
              call mxm(dxm12,nx2,x(1,e),nx1,ta1,nyz1)
              call mxm(ta1,nx2,iytm12,ny1,dx(1,e),ny2)
              call col2(dx(1,e),rm2(1,e),nxyz2)
              call mxm(ixm12,nx2,x(1,e),nx1,ta1,nyz1)
              call mxm(ta1,nx2,dytm12,ny1,ta3,ny2)
              call addcol3(dx(1,e),ta3,sm2(1,e),nxyz2)
          endif
#endif
      else  ! 3D

      !           if (ifsplit) then
      
      !             call mxm  (dxm12,nx2,x(1,e),nx1,dx(1,e),nyz1)
      !             call col2 (dx(1,e),rm2(1,e),nxyz2)
      !             i1=1
      !             i2=1
      !             do iz=1,nz1
      !                call mxm (x(1,e),nx2,dytm12,ny1,ta1(i2),ny2)
      !                i1=i1+n1
      !                i2=i2+n2
      !             enddo
      !             call addcol3 (dx(1,e),ta1,sm2(1,e),nxyz2)
      !             call mxm (x(1,e),nxy2,dztm12,nz1,ta1,nz2)
      !             call addcol3 (dx(1,e),ta1,tm2(1,e),nxyz2)

      !           else ! PN - PN-2

          call mxm(dxm12,nx2,x(:,:,:,e),nx1,ta1,nyz1)
          i1=1
          i2=1
          do iz=1,nz1
              call mxm(ta1(i1),nx2,iytm12,ny1,ta2(i2),ny2)
              i1=i1+n1
              i2=i2+n2
          enddo
          call mxm(ta2,nxy2,iztm12,nz1,dx(:,:,:,e),nz2)
          dx(:,:,:,e) = dx(:,:,:,e) * rm2(:,:,:,e)

          call mxm(ixm12,nx2,x(:,:,:,e),nx1,ta3,nyz1) ! reuse ta3 below
          i1=1
          i2=1
          do iz=1,nz1
              call mxm(ta3(i1),nx2,dytm12,ny1,ta2(i2),ny2)
              i1=i1+n1
              i2=i2+n2
          enddo
          call mxm(ta2,nxy2,iztm12,nz1,ta1,nz2)
          dx(:,:,:,e) = dx(:,:,:,e) + reshape(ta1,(/lx2,ly2,lz2/)) * sm2(:,:,:,e)

      !            call mxm (ixm12,nx2,x(1,e),nx1,ta1,nyz1) ! reuse ta3 from above
          i1=1
          i2=1
          do iz=1,nz1
              call mxm(ta3(i1),nx2,iytm12,ny1,ta2(i2),ny2)
              i1=i1+n1
              i2=i2+n2
          enddo
          call mxm(ta2,nxy2,dztm12,nz1,ta3,nz2)
          dx(:,:,:,e) = dx(:,:,:,e) + reshape(ta3,(/lx2,ly2,lz2/)) * tm2(:,:,:,e)
      !           endif
      endif
  
  !        Collocate with the weights on the pressure mesh


      if(ifsplit) then
          dx(:,:,:,e) = dx(:,:,:,e) * bm1(:,:,:,e)
          dx(:,:,:,e) = dx(:,:,:,e) / jacm1(:,:,:,e)
      else
          write(*,*) "Whoops! multd"
#if 0
          if ( .NOT. ifaxis) call col2(dx(1,e),w3m2,nxyz2)
          if (ifaxis) then
              if (ifrzer(e)) then
                  call col2(dx(1,e),bm2(1,1,1,e),nxyz2)
                  dx(:,e) = dx(:,e) / reshape(jacm2(:,:,:,e), (/nxyz2/))
              else
                  call col2(dx(1,e),w3m2,nxyz2)
                  call col2(dx(1,e),ym2(1,1,1,e),nxyz2)
              endif
          endif
#endif
      endif

  !        If axisymmetric, add an extra diagonal term in the radial
  !        direction (ISD=2).
  !        NOTE: NZ1=NZ2=1

      if(ifsplit) then

          if (ifaxis .AND. (isd == 2) .AND. iflg == 1) then
              write(*,*) "Whoops! multd"
#if 0
              call copy(ta3,x(1,e),nxyz1)
              if (ifrzer(e)) then
                  call rzero(ta3, nx1)
                  call mxm(x(1,e),nx1,datm1,ny1,duax,1)
                  call copy(ta3,duax,nx1)
              endif
              call col2(ta3,baxm1(1,1,1,e),nxyz1)
              call add2(dx(1,e),ta3,nxyz2)
#endif
          endif

      else
          write(*,*) "Whoops! multd"
#if 0
          if (ifaxis .AND. (isd == 2)) then
              call mxm(ixm12,nx2,x(1,e),nx1,ta1,ny1)
              call mxm(ta1,nx2,iytm12,ny1,ta2,ny2)
              call col3(ta3,bm2(1,1,1,e),ta2,nxyz2)
              call invcol2(ta3,ym2(1,1,1,e),nxyz2)
              call add2(dx(1,e),ta3,nxyz2)
          endif
#endif
      endif
  
  enddo

#ifndef NOTIMER
  tmltd=tmltd+(dnekclock()-etime1)
#endif
  return
end subroutine multd

!----------------------------------------------------------------------
!----------------------------------------------------------------------
subroutine ophx (out1,out2,out3,inp1,inp2,inp3,h1,h2)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, ndim
  use input, only : ifstrs
  implicit none

  REAL(DP) :: OUT1 (lx1,ly1,lz1,1)
  REAL(DP) :: OUT2 (lx1,ly1,lz1,1)
  REAL(DP) :: OUT3 (lx1,ly1,lz1,1)
  REAL(DP) :: INP1 (lx1,ly1,lz1,1)
  REAL(DP) :: INP2 (lx1,ly1,lz1,1)
  REAL(DP) :: INP3 (lx1,ly1,lz1,1)
  REAL(DP) :: H1   (lx1,ly1,lz1,1)
  REAL(DP) :: H2   (lx1,ly1,lz1,1)

  integer :: imesh
  imesh = 1

  IF (ifstrs) THEN
#if 0
      matmod = 0
      CALL axhmsf(out1,out2,out3,inp1,inp2,inp3,h1,h2,matmod)
#endif
  ELSE
      CALL axhelm(out1,inp1,h1,h2,imesh,1)
      CALL axhelm(out2,inp2,h1,h2,imesh,2)
      IF (ndim == 3) &
      CALL axhelm(out3,inp3,h1,h2,imesh,3)
  ENDIF

  return
end subroutine ophx

!--------------------------------------------------------------
!--------------------------------------------------------------
subroutine dudxyz (du,u,rm1,sm1,tm1,jm1,imsh,isd)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1
  use size_m, only : nx1, ny1, nz1, nelv, nelt, ndim
  use dxyz, only : dxm1, dytm1, dztm1
  use geom, only : jacmi
  implicit none

  integer :: imsh, isd
  REAL(DP) ::  DU  (lx1,ly1,lz1,*)
  REAL(DP) ::  U   (lx1,ly1,lz1,*)
  REAL(DP) ::  RM1 (lx1,ly1,lz1,*)
  REAL(DP) ::  SM1 (lx1,ly1,lz1,*)
  REAL(DP) ::  TM1 (lx1,ly1,lz1,*)
  REAL(DP) ::  JM1 (lx1,ly1,lz1,*)

  REAL(DP) ::  DRST(lx1,ly1,lz1)
  integer :: nel, nxy1, nyz1, nxyz1, ntot, iel, iz

  nel = -1
  IF (imsh == 1) nel = nelv
  IF (imsh == 2) nel = nelt
  nxy1  = nx1*ny1
  nyz1  = ny1*nz1
  nxyz1 = nx1*ny1*nz1
  ntot  = nxyz1*nel

  DO 1000 iel=1,nel
  
!max      IF (IFAXIS) CALL SETAXDY (IFRZER(IEL) )
  
      IF (ndim == 2) THEN
          CALL mxm(dxm1,nx1,u(1,1,1,iel),nx1,du(1,1,1,iel),nyz1)
          du(:,:,:,iel) = du(:,:,:,iel) * rm1(:,:,:,iel)
          CALL mxm(u(1,1,1,iel),nx1,dytm1,ny1,drst,ny1)
          du(:,:,:,iel) = du(:,:,:,iel) + drst * sm1(:,:,:,iel)
      ELSE
          CALL mxm(dxm1,nx1,u(1,1,1,iel),nx1,du(1,1,1,iel),nyz1)
          du(:,:,:,iel) = du(:,:,:,iel) * rm1(:,:,:,iel)
          DO 20 iz=1,nz1
              CALL mxm(u(1,1,iz,iel),nx1,dytm1,ny1,drst(1,1,iz),ny1)
          20 END DO
          du(:,:,:,iel) = du(:,:,:,iel) + drst * sm1(:,:,:,iel)
          CALL mxm(u(1,1,1,iel),nxy1,dztm1,nz1,drst,nz1)
          du(:,:,:,iel) = du(:,:,:,iel) + drst * tm1(:,:,:,iel)
      ENDIF
  
  1000 END DO

!    CALL INVCOL2 (DU,JM1,NTOT)
  du(:,:,:,1:nel) = du(:,:,:,1:nel) * jacmi

  return
end subroutine dudxyz

!---------------------------------------------------------------------
!----------------------------------------------------------------------
subroutine makef
  use kinds, only : dp
  use input, only : ifnav, ifchar, iftran
  use ctimer, only : tmakef, nmakef, dnekclock
  implicit none
  real(DP) :: etime

  nmakef = nmakef + 1
  etime = dnekclock()
  call makeuf
!  if (ifnatc)                               call natconv
!  if (ifexplvis .AND. ifsplit)                call explstrs
  etime = etime - dnekclock()
  if (ifnav .AND. ( .NOT. ifchar))              call advab
  etime = etime + dnekclock()
!  if (ifmvbd)                               call admeshv
  if (iftran) then
    call makeabf
  endif
  if ((iftran .AND. .NOT. ifchar) .OR. &
  (iftran .AND. .NOT. ifnav .AND. ifchar))   call makebdf
!max    if (ifnav .AND. ifchar .AND. ( .NOT. ifmvbd))   call advchar
#if 0
  if (ifmodel)                              call twallsh
#endif
  tmakef = tmakef + (dnekclock() - etime)
  return
end subroutine makef

!---------------------------------------------------------------------
!----------------------------------------------------------------------
subroutine makeuf
  use size_m, only : nelv
  use geom, only : bm1
  use soln, only : bfx, bfy, bfz
  use tstep, only : time, dt
  implicit none

  integer :: i
  time = time-dt
  CALL nekuf(bfx,bfy,bfz)
  do i = 1, nelv
    bfx(:,:,:,i) = bfx(:,:,:,i) * bm1(:,:,:,i)
    bfy(:,:,:,i) = bfy(:,:,:,i) * bm1(:,:,:,i)
    bfz(:,:,:,i) = bfz(:,:,:,i) * bm1(:,:,:,i)
  enddo
  time = time+dt

  return
end subroutine makeuf

subroutine nekuf (f1,f2,f3)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lelv, nx1, ny1, nz1, nelv
  use nekuse, only : ffx, ffy, ffz
  use parallel, only : lglel
  implicit none

  REAL(DP) :: F1 (lx1,ly1,lz1,lelv)
  REAL(DP) :: F2 (lx1,ly1,lz1,lelv)
  REAL(DP) :: F3 (lx1,ly1,lz1,lelv)

  integer :: i, j, k, ielg, iel
  f1 = 0._dp; f2 = 0._dp; f3=0._dp 
  DO iel=1,nelv
      ielg = lglel(iel)
      DO k=1,nz1
          DO j=1,ny1
              DO i=1,nx1
                  CALL nekasgn(i,j,k,iel)
                  CALL userf(i,j,k,ielg)
                  f1(i,j,k,iel) = ffx
                  f2(i,j,k,iel) = ffy
                  f3(i,j,k,iel) = ffz
              enddo
          enddo
      enddo
  END DO

  return
end subroutine nekuf

!---------------------------------------------------------------
!---------------------------------------------------------------
subroutine advab()
  use kinds, only : dp
  use size_m, only : nx1, ny1, nz1, nelv, ndim
  use geom, only : bm1
  use soln, only : vx, vy, vz, bfx, bfy, bfz
  implicit none

  real(DP), allocatable:: TA (:,:,:,:)

  integer :: ntot1

  allocate(ta(nx1,ny1,nz1,nelv)) 
  ntot1 = nx1*ny1*nz1*nelv
  CALL convop(ta,vx)
  bfx = bfx - bm1*ta
  CALL convop(ta,vy)
  bfy = bfy - bm1*ta
  IF (ndim /= 2) THEN
      CALL convop(ta,vz)
      bfz = bfz - bm1*ta
  ENDIF

  return
end subroutine advab

!-----------------------------------------------------------------------
subroutine makebdf()
  use kinds, only : dp
  use size_m, only : nx1, ny1, nz1, nelv
  use geom, only : ifgeom
  use geom, only : bm1, bm1lag
  use soln, only : vx, vy, vz, vxlag, vylag, vzlag, bfx, bfy, bfz
  use soln, only : vtrans
  use tstep, only : dt, ifield, bd, nbd
  implicit none

  real(DP), allocatable :: H2 (:,:,:)

  integer :: ilag, i

  allocate(h2(nx1,ny1,nz1))
 
  do i = 1, nelv 
    h2 = (1./dt) * vtrans(:,:,:,i,ifield) 
    DO ilag=2,nbd
      IF (ifgeom) THEN
        bfx(:,:,:,i) = bfx(:,:,:,i) + h2 * (vxlag(:,:,:,i,ilag-1) * bm1lag(:,:,:,i,ilag-1) * bd(ilag+1))
        bfy(:,:,:,i) = bfy(:,:,:,i) + h2 * (vylag(:,:,:,i,ilag-1) * bm1lag(:,:,:,i,ilag-1) * bd(ilag+1))
        bfz(:,:,:,i) = bfz(:,:,:,i) + h2 * (vzlag(:,:,:,i,ilag-1) * bm1lag(:,:,:,i,ilag-1) * bd(ilag+1))
      ELSE
        bfx(:,:,:,i) = bfx(:,:,:,i) + h2 * (vxlag(:,:,:,i,ilag-1) * bm1(:,:,:,i)    * bd(ilag+1))
        bfy(:,:,:,i) = bfy(:,:,:,i) + h2 * (vylag(:,:,:,i,ilag-1) * bm1(:,:,:,i)    * bd(ilag+1))
        bfz(:,:,:,i) = bfz(:,:,:,i) + h2 * (vzlag(:,:,:,i,ilag-1) * bm1(:,:,:,i)    * bd(ilag+1))
      ENDIF
    END DO
    bfx(:,:,:,i) = bfx(:,:,:,i) + h2*vx(:,:,:,i)*bm1(:,:,:,i)*bd(2)
    bfy(:,:,:,i) = bfy(:,:,:,i) + h2*vy(:,:,:,i)*bm1(:,:,:,i)*bd(2)
    bfz(:,:,:,i) = bfz(:,:,:,i) + h2*vz(:,:,:,i)*bm1(:,:,:,i)*bd(2)
  enddo

  return
end subroutine makebdf

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine makeabf
  use kinds, only : dp
  use size_m, only : nx1, ny1, nz1, nelv, ndim
  use soln, only : abx1, aby1, abz1, abx2, aby2, abz2, bfx, bfy, bfz, vtrans
  use tstep, only : ab
  implicit none

  real(DP), allocatable :: TA (:,:,:,:) 

  integer :: ntot1, iel
  real(DP) :: ab0, ab1, ab2

  allocate(ta(nx1,ny1,nz1,nelv))

  ntot1 = nx1*ny1*nz1*nelv

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

#if 0
! 11*ntot mops
! 6*ntot flops
  ta = ab1 * abx1 + ab2 * abx2
  CALL copy(abx2,abx1,ntot1)
  CALL copy(abx1,bfx,ntot1)
  bfx = (ab0 * bfx + ta) * vtrans(:,:,:,:,1) ! multiply by density

  ta = ab1 * aby1 + ab2 * aby2
  CALL copy(aby2,aby1,ntot1)
  CALL copy(aby1,bfy,ntot1)
  bfy = (ab0 * bfy + ta) * vtrans(:,:,:,:,1) ! multiply by density

  IF (ndim == 3) THEN
    ta = ab1 * abz1 + ab2 * abz2
    CALL copy(abz2,abz1,ntot1)
    CALL copy(abz1,bfz,ntot1)
    bfz = (ab0 * bfz + ta) * vtrans(:,:,:,:,1) ! multiply by density
  ENDIF
#else
! 7*ntot mops
! 6*ntot flops
  do iel = 1, nelv
    ta(:,:,:,1) = ab1 * abx1(:,:,:,iel) + ab2 * abx2(:,:,:,iel)
    abx2(:,:,:,iel) = abx1(:,:,:,iel)
    abx1(:,:,:,iel) = bfx(:,:,:,iel)
    bfx(:,:,:,iel) = (ab0*bfx(:,:,:,iel) + ta(:,:,:,1)) * vtrans(:,:,:,iel,1)
  enddo

  do iel = 1, nelv
    ta(:,:,:,1) = ab1 * aby1(:,:,:,iel) + ab2 * aby2(:,:,:,iel)
    aby2(:,:,:,iel) = aby1(:,:,:,iel)
    aby1(:,:,:,iel) = bfy(:,:,:,iel)
    bfy(:,:,:,iel) = (ab0*bfy(:,:,:,iel) + ta(:,:,:,1)) * vtrans(:,:,:,iel,1)
  enddo

  do iel = 1, nelv
    ta(:,:,:,1) = ab1 * abz1(:,:,:,iel) + ab2 * abz2(:,:,:,iel)
    abz2(:,:,:,iel) = abz1(:,:,:,iel)
    abz1(:,:,:,iel) = bfz(:,:,:,iel)
    bfz(:,:,:,iel) = (ab0*bfz(:,:,:,iel) + ta(:,:,:,1)) * vtrans(:,:,:,iel,1)
  enddo
#endif

  return
end subroutine makeabf

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine setabbd (ab,dtlag,nab,nbd)
  use kinds, only : dp
  implicit none

  REAL(DP), intent(out) :: AB(nab)    !>
  real(DP), intent(in)  :: DTLAG(nbd) !>
  integer,  intent(in)  :: nab        !>
  integer,  intent(in)  :: nbd        !>

  real(DP) :: dt0, dt1, dt2, dta, dts, dtb, dtc, dtd, dte

  IF ( nab == 1 ) THEN
  
      ab(1) = 1.0
  
  ELSEIF ( nab == 2 ) THEN
      dt0 = dtlag(1)
      dt1 = dtlag(2)
  
      dta =  dt0/dt1
  
      IF ( nbd == 1 ) THEN
      
          ab(2) = -0.5*dta
          ab(1) =  1.0 - ab(2)
      
      ELSEIF ( nbd == 2 ) THEN
      
          ab(2) = -dta
          ab(1) =  1.0 - ab(2)
      
      ENDIF
  
  ELSEIF ( nab == 3 ) THEN
      dt0 = dtlag(1)
      dt1 = dtlag(2)  
      dt2 = dtlag(3)

      dts =  dt1 + dt2
      dta =  dt0 / dt1
      dtb =  dt1 / dt2
      dtc =  dt0 / dt2
      dtd =  dts / dt1
      dte =  dt0 / dts
  
      IF ( nbd == 1 ) THEN
      
          ab(3) =  dte*( 0.5*dtb + dtc/3. )
          ab(2) = -0.5*dta - ab(3)*dtd
          ab(1) =  1.0 - ab(2) - ab(3)
      
      ELSEIF ( nbd == 2 ) THEN
      
          ab(3) =  2./3.*dtc*(1./dtd + dte)
          ab(2) = -dta - ab(3)*dtd
          ab(1) =  1.0 - ab(2) - ab(3)
      
      ELSEIF ( nbd == 3 ) THEN
      
          ab(3) =  dte*(dtb + dtc)
          ab(2) = -dta*(1.0 + dtb + dtc)
          ab(1) =  1.0 - ab(2) - ab(3)
      
      ENDIF
  
  ENDIF

  return
end subroutine setabbd

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine setbd (bd,dtbd,nbd)
  use kinds, only : dp
  implicit none

  REAL(dp), intent(out) :: BD(*)   !>
  real(dp), intent(in)  :: DTBD(*) !>
  integer , intent(in)  :: nbd     !>

  integer, PARAMETER :: NDIM = 10
  REAL(DP) :: BDMAT(ndim,ndim),BDRHS(ndim), BDF
  INTEGER :: IR(ndim),IC(ndim)
  integer :: nsys, i, ibd

  bd(1:ndim) = 0._dp; bdf = -1
  ! BDF(1) is trivial
  IF (nbd == 1) THEN
      bd(1) = 1.
      bdf   = 1.
  ! BDF(>1) computed using a linear system
  ELSEIF (nbd >= 2) THEN
      nsys = nbd+1
      CALL bdsys(bdmat,bdrhs,dtbd,nbd,ndim)
      CALL lu(bdmat,nsys,ndim,ir,ic)
      CALL solve(bdrhs,bdmat,1,nsys,ndim,ir,ic)
      DO 30 i=1,nbd
          bd(i) = bdrhs(i)
      30 END DO
      bdf = bdrhs(nbd+1)
  ENDIF

  !   Normalize
  DO ibd=nbd,1,-1
      bd(ibd+1) = bd(ibd)
  END DO
  bd(1) = 1.
  DO ibd=1,nbd+1
      bd(ibd) = bd(ibd)/bdf
  END DO
!   write(6,1) (bd(k),k=1,nbd+1)
! 1 format('bd:',1p8e13.5)

  return
end subroutine setbd

subroutine bdsys (a,b,dt,nbd,ndim)
  use kinds, only : dp
  implicit none

  integer :: nbd, ndim
  REAL(DP) :: A(ndim,9),B(9),DT(9)

  integer :: n, j, k, i
  real(DP) :: sumdt
  a = 0._dp
  n = nbd+1
  DO j=1,nbd
      a(1,j) = 1.
  END DO
  a(1,nbd+1) = 0.
  b(1) = 1.
  DO j=1,nbd
      sumdt = 0.
      DO k=1,j
          sumdt = sumdt+dt(k)
      END DO
      a(2,j) = sumdt
  END DO
  a(2,nbd+1) = -dt(1)
  b(2) = 0.
  DO i=3,nbd+1
      DO j=1,nbd
          sumdt = 0.
          DO k=1,j
              sumdt = sumdt+dt(k)
          END DO
          a(i,j) = sumdt**(i-1)
      END DO
      a(i,nbd+1) = 0.
      b(i) = 0.
  END DO
  return
end subroutine bdsys

!-------------------------------------------------------------------
!-------------------------------------------------------------------
subroutine tauinit (tau,ilag)
  use kinds, only : dp
  use tstep, only : nbd, dtlag
  implicit none

  real(DP) :: tau
  integer :: ilag
  integer :: i

  tau   = 0.
  DO 10 i=nbd,ilag+1,-1
      tau = tau+dtlag(i)
  10 END DO
  return
end subroutine tauinit

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine lagvel
  use size_m, only : nx1, ny1, nz1, nelv, ndim
  use soln, only : vxlag, vylag, vzlag, vx, vy, vz
  implicit none

  integer :: ntot1, ilag
  ntot1 = nx1*ny1*nz1*nelv

!    DO 100 ILAG=NBDINP-1,2,-1
  DO 100 ilag=3-1,2,-1
      CALL copy(vxlag(1,1,1,1,ilag),vxlag(1,1,1,1,ilag-1),ntot1)
      CALL copy(vylag(1,1,1,1,ilag),vylag(1,1,1,1,ilag-1),ntot1)
      IF (ndim == 3) &
      CALL copy(vzlag(1,1,1,1,ilag),vzlag(1,1,1,1,ilag-1),ntot1)
  100 END DO

  CALL opcopy(vxlag,vylag,vzlag,vx,vy,vz)

  return
end subroutine lagvel

!----------------------------------------------------------------------
!----------------------------------------------------------------------
subroutine setordbd
  use tstep, only : nbdinp, nbd, istep
  implicit none

!   IF (IFSPLIT .OR. NBDINP.EQ.0) THEN     undid hardwire, 3/6/92 pff
  IF ( nbdinp < 1) THEN
      nbd = 1
  ELSE
      IF ((istep == 0) .OR. (istep == 1))        nbd = 1
      IF ((istep > 1) .AND. (istep <= nbdinp))  nbd = istep
      IF (istep > nbdinp)                     nbd = nbdinp
  ENDIF

  return
end subroutine setordbd

!---------------------------------------------------------------
!---------------------------------------------------------------
subroutine normsc (h1,semi,l2,linf,x,imesh)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lelt
  use size_m, only : nx1, ny1, nz1, nelv, nelt, ndim
  use geom, only : bm1, voltm1, volvm1
  use ctimer, only : tnmsc, nnmsc, dnekclock
  implicit none

  real(DP) :: h1, semi, l2, linf
  real(DP), intent(in) :: X  (lx1,ly1,lz1,*)
  integer :: imesh

  real(DP), allocatable, dimension(:,:,:,:) :: y, ta1, ta2
  REAL(DP) :: LENGTH, vol
  integer :: nel, nxyz1, ntot1
  real(DP), external :: glamax, glsum
  real(DP) :: etime

  nnmsc = nnmsc + 1
  etime = dnekclock()

  allocate(y(lx1,ly1,lz1,lelt), ta1(lx1,ly1,lz1,lelt), ta2(lx1,ly1,lz1,lelt))

  IF (imesh == 1) THEN
      nel = nelv
      vol = volvm1
  ELSEIF (imesh == 2) THEN
      nel = nelt
      vol = voltm1
  else
    nel = -1
    vol = -1.
    write(*,*) "IMESH \notin {1,2}"
  ENDIF

  length = vol**(1./(ndim))
  nxyz1  = nx1*ny1*nz1
  ntot1  = nxyz1*nel

  h1     = 0.
  semi   = 0.
  l2     = 0.
  linf   = 0.

  linf = glamax(x,ntot1)

  ta1 = x(:,:,:,1:nel)*x(:,:,:,1:nel) * bm1
  l2   = glsum(ta1,ntot1)
  IF (l2 < 0.0) l2 = 0.

  ta1 = 1._dp; ta2 = 0._dp
  etime = etime - dnekclock()
  CALL axhelm(y,x,ta1,ta2,imesh,1)
  etime = etime + dnekclock()
  ta1 = y * x(:,:,:,1:nel)
  semi = glsum(ta1,ntot1)
  IF (semi < 0.0) semi = 0.

  h1   = sqrt((semi*length**2+l2)/vol)
  semi = sqrt(semi/vol)
  l2   = sqrt(l2/vol)
  IF (h1 < 0.) h1 = 0.
  tnmsc = tnmsc + (dnekclock() - etime)

  return
end subroutine normsc

!---------------------------------------------------------------
! defined on mesh 1 (velocity mesh only !)
! The error norms are normalized with respect to the volume
! (except for Linf).
!---------------------------------------------------------------
subroutine normvc (h1,semi,l2,linf,x1,x2,x3)
  use kinds, only : dp
  use size_m, only : nx1, ny1, nz1, nelv, ndim
  use size_m, only : lx1, ly1, lz1, lelt
  use geom, only : volvm1, bm1
  use ctimer, only : tnmvc, nnmvc, dnekclock
  implicit none

  REAL(DP) :: H1,SEMI,L2,LINF
  REAL(DP) :: X1 (lx1,ly1,lz1,lelt)
  REAL(DP) :: X2 (lx1,ly1,lz1,lelt)
  REAL(DP) :: X3 (lx1,ly1,lz1,lelt)

  real(DP), allocatable, dimension(:,:,:,:) :: Y1, Y2, Y3, TA1, TA2
  REAL(DP) :: LENGTH
  integer :: imesh, nel, nxyz1, ntot1
  real(DP) :: vol
  real(DP), external :: glamax, glsum
  real(DP) :: etime

  nnmvc = nnmvc + 1
  etime = dnekclock()

  imesh  = 1
  nel    = nelv
  vol    = volvm1
  length = vol**(1./(ndim))
  nxyz1  = nx1*ny1*nz1
  ntot1  = nxyz1*nel

  h1     = 0.
  semi   = 0.
  l2     = 0.
  linf   = 0.

  allocate(y1(lx1,ly1,lz1,lelt) &
  ,y2(lx1,ly1,lz1,lelt) &
  ,y3(lx1,ly1,lz1,lelt) &
  ,ta1(lx1,ly1,lz1,lelt) &
  ,ta2(lx1,ly1,lz1,lelt) )


  IF (ndim == 3) THEN
    ta1 = x1*x1 + x2*x2 + x3*x3
  else
    ta1 = x1*x1 + x2*x2
  ENDIF
  linf = glamax(ta1,ntot1)
  linf = sqrt( linf )

  ta1 = ta1 * bm1
  l2 = glsum(ta1,ntot1)
  IF (l2 < 0.0) l2 = 0.

  ta1 = 1._dp; ta2 = 0._dp
  etime = etime - dnekclock()
  CALL ophx(y1,y2,y3,x1,x2,x3,ta1,ta2)
  etime = etime + dnekclock()
  IF (ndim == 3) THEN
    ta1 = y1*x1 + y2*x2 + y3*x3
  else
    ta1 = y1*x1 + y2*x2
  ENDIF

  semi = glsum(ta1,ntot1)
  IF (semi < 0.0) semi = 0.

  h1   = sqrt((semi*length**2+l2)/vol)
  semi = sqrt(semi/vol)
  l2   = sqrt(l2  /vol)
  IF (h1 < 0.) h1 = 0.

  tnmvc = tnmvc + (dnekclock() - etime)

  return
end subroutine normvc

subroutine opcolv (a1,a2,a3,c)
  use kinds, only : dp
  use size_m, only : nx1, ny1, nz1, nelv, ndim
  use opctr, only : isclld, nrout, myrout, rname, dct, ncall, dcount
  implicit none
  REAL(DP) :: A1(1),A2(1),A3(1),C(1)
  integer :: ntot1, i, isbcnt

  ntot1=nx1*ny1*nz1*nelv

#ifndef NOTIMER
  if (isclld == 0) then
      isclld=1
      nrout=nrout+1
      myrout=nrout
      rname(myrout) = 'opcolv'
  endif

  isbcnt = ntot1*ndim
  dct(myrout) = dct(myrout) + (isbcnt)
  ncall(myrout) = ncall(myrout) + 1
  dcount      =      dcount + (isbcnt)
#endif

  IF (ndim == 3) THEN
      DO 100 i=1,ntot1
          a1(i)=a1(i)*c(i)
          a2(i)=a2(i)*c(i)
          a3(i)=a3(i)*c(i)
      100 END DO
  ELSE
      DO 200 i=1,ntot1
          a1(i)=a1(i)*c(i)
          a2(i)=a2(i)*c(i)
      200 END DO
  ENDIF
  return
end subroutine opcolv

subroutine opcopy (a1,a2,a3,b1,b2,b3)
  use kinds, only : dp
  use size_m
  implicit none
  REAL(DP) :: A1(1),A2(1),A3(1),B1(1),B2(1),B3(1)
  integer :: ntot1
  ntot1=nx1*ny1*nz1*nelv
  CALL copy(a1,b1,ntot1)
  CALL copy(a2,b2,ntot1)
  IF(ndim == 3)CALL copy(a3,b3,ntot1)
  return
end subroutine opcopy

!-----------------------------------------------------------------------
subroutine opdssum (a,b,c)! NOTE: opdssum works on FLUID/MHD arrays only!
  use kinds, only : dp
  use input, only : ifcyclic
  implicit none

  real(DP) :: a(1),b(1),c(1)

  if (ifcyclic) then
    write(*,*) "Oops: ifcyclic"
#if 0
      call rotate_cyc(a,b,c,1)
      call vec_dssum(a,b,c,nx1,ny1,nz1)
      call rotate_cyc(a,b,c,0)
#endif
  else
      call vec_dssum(a,b,c)
  endif

  return
end subroutine opdssum

!-----------------------------------------------------------------------
subroutine opdsop (a,b,c,op)! opdsop works on FLUID/MHD arrays only!
  use kinds, only : dp
  use input, only : ifcyclic
  implicit none

  real(DP) :: a(1),b(1),c(1)
  character(3) :: op

  if (ifcyclic) then
      if (op == '*  ' .OR. op == 'mul' .OR. op == 'MUL') then
          call vec_dsop(a,b,c,op)
      else
          write(*,*) "Oops: op"
#if 0
          call rotate_cyc(a,b,c,1)
          call vec_dsop(a,b,c,nx1,ny1,nz1,op)
          call rotate_cyc(a,b,c,0)
#endif
      endif
  else
      call vec_dsop(a,b,c,op)
  endif

  return
end subroutine opdsop

!-----------------------------------------------------------------------
subroutine transpose(a,lda,b,ldb)
  use kinds, only : dp
  implicit none
  integer :: lda, ldb
  real(DP) :: a(lda,*),b(ldb,*)
  integer :: i, j

  do j=1,ldb
      do i=1,lda
          a(i,j) = b(j,i)
      enddo
  enddo
  return
end subroutine transpose

!-----------------------------------------------------------------------
! using the skew-symmetric formulation.
! The field variable FI is defined on mesh M1 (GLL) and
! the velocity field is assumed given.
subroutine convop(conv,fi)
  use kinds, only : dp
  use ctimer, only : icalld, tadvc, nadvc, etime1, dnekclock
  use size_m, only : lx1, ly1, lz1
  use size_m, only : nx1, ny1, nz1, nelv
  use dealias, only : vxd, vyd, vzd
  use input, only : param, ifpert
  use geom, only : bm1
  use soln, only : vx, vy, vz
  use tstep, only : nelfld, ifield
  implicit none

!     Arrays in parameter list
  REAL(DP) :: CONV (lx1,ly1,lz1,*)
  REAL(DP) :: FI   (lx1,ly1,lz1,*)

  integer :: nxyz1, ntot1, ntotz

#ifndef NOTIMER
  nadvc=nadvc + 1
  etime1=dnekclock()
#endif

  nxyz1 = nx1*ny1*nz1
  ntot1 = nx1*ny1*nz1*nelv
  ntotz = nx1*ny1*nz1*nelfld(ifield)

  conv(:,:,:,1:nelfld(ifield)) = 0._dp

  if (param(86) /= 0.0) then  ! skew-symmetric form
!max        call convopo(conv,fi)
      goto 100
  endif

!     write(6,*) istep,param(99),' CONVOP',ifpert
!     ip99 = param(99)
!     if (istep.gt.5) call exitti(' CONVOP dbg: $',ip99)

  if (param(99) == 2 .OR. param(99) == 3) then
!max        call conv1d(conv,fi)  !    use dealiased form
  elseif (param(99) == 4) then
      if (ifpert) then
          etime1 = etime1 - dnekclock()
          call convect_new(conv,fi, .false. ,vx,vy,vz, .false. )
          etime1 = etime1 + dnekclock()
      else
          etime1 = etime1 - dnekclock()
          call convect_new(conv,fi, .false. ,vxd,vyd,vzd, .true. )
          etime1 = etime1 + dnekclock()
      endif
      conv(:,:,:,1:nelv) = conv(:,:,:,1:nelv) / bm1 ! local mass inverse
  elseif (param(99) == 5) then
!max        call convect_cons(conv,fi, .FALSE. ,vx,vy,vz, .FALSE. )
      conv(:,:,:,1:nelv) = conv(:,:,:,1:nelv) / bm1  ! local mass inverse
  else
!max        call conv1 (conv,fi)  !    use the convective form
  endif

  100 continue

#ifndef NOTIMER
  tadvc=tadvc+(dnekclock()-etime1)
#endif

  return
end subroutine convop

!-----------------------------------------------------------------------
!---------------------------------------------------------------------
subroutine opdiv(outfld,inpx,inpy,inpz)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lx2, ly2, lz2, lelv
  use size_m, only : nx2, ny2, nz2, nelv, ndim, nx1, ny1, nz1
  use geom, only : rxm2, rym2, rzm2, sxm2, sym2, szm2, txm2, tym2, tzm2
  use mesh, only : if_ortho
  use geom, only : bm1, jacmi
  use dxyz, only : dxm12, dytm12, dztm12
  use ixyz, only : ixm12, iytm12, iztm12

  implicit none

  real(DP) :: outfld (lx2,ly2,lz2,nelv)
  real(DP) :: inpx   (lx1,ly1,lz1,nelv)
  real(DP) :: inpy   (lx1,ly1,lz1,nelv)
  real(DP) :: inpz   (lx1,ly1,lz1,nelv)
  
  real(DP), allocatable :: work (:,:,:,:)

  real(DP) ::  ta1 (lx1*ly1*lz1) &
  ,             ta2 (lx1*ly1*lz1) &
  ,             ta3 (lx1*ly1*lz1)
  integer :: iflg, ntot2, i1, i2, e, iz, n1, n2, nxy2, nyz1

  allocate(work(lx2,ly2,lz2,lelv))

  iflg = 1

  ntot2 = nx2*ny2*nz2*nelv
  if (if_ortho) then
    nxy2  = nx2*ny2
    n1    = nx2*ny1
    n2    = nx2*ny2
    nyz1  = ny1*nz1
    do e=1,nelv
      call mxm(dxm12,nx2,inpx(:,:,:,e),nx1,ta1,nyz1)
      i1=1
      i2=1
      do iz=1,nz1
          call mxm(ta1(i1),nx2,iytm12,ny1,ta2(i2),ny2)
          i1=i1+n1
          i2=i2+n2
      enddo
      call mxm(ta2,nxy2,iztm12,nz1,outfld(:,:,:,e),nz2)
      outfld(:,:,:,e) = outfld(:,:,:,e) * rxm2(:,:,:,e)
 
      call mxm(ixm12,nx2,inpy(:,:,:,e),nx1,ta3,nyz1) ! reuse ta3 below
      i1=1
      i2=1
      do iz=1,nz1
          call mxm(ta3(i1),nx2,dytm12,ny1,ta2(i2),ny2)
          i1=i1+n1
          i2=i2+n2
      enddo
      call mxm(ta2,nxy2,iztm12,nz1,ta1,nz2)
      outfld(:,:,:,e) = outfld(:,:,:,e) + reshape(ta1,(/lx2,ly2,lz2/)) * sym2(:,:,:,e)
 
      call mxm(ixm12,nx2,inpz(:,:,:,e),nx1,ta1,nyz1) 
      i1=1
      i2=1
      do iz=1,nz1
          call mxm(ta3(i1),nx2,iytm12,ny1,ta2(i2),ny2)
          i1=i1+n1
          i2=i2+n2
      enddo
      call mxm(ta2,nxy2,dztm12,nz1,ta3,nz2)
      outfld(:,:,:,e) = outfld(:,:,:,e) + reshape(ta3,(/lx2,ly2,lz2/)) * tzm2(:,:,:,e)
    enddo
    outfld = outfld * bm1 * jacmi

  else
    call multd(work,inpx,rxm2,sxm2,txm2,1,iflg)
    call copy(outfld,work,ntot2)
    call multd(work,inpy,rym2,sym2,tym2,2,iflg)
    outfld = outfld + work
    if (ndim == 3) then
        call multd(work,inpz,rzm2,szm2,tzm2,3,iflg)
        outfld = outfld + work
    endif
  endif

  return
end subroutine opdiv

!-----------------------------------------------------------------------
subroutine wgradm1(ux,uy,uz,u,nel) ! weak form of grad
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, nx1
  use dxyz, only : dxm1, dxtm1
  use geom, only : rxm1, sxm1, txm1, rym1, sym1, tym1, rzm1, szm1, tzm1
  use input, only : if3d
  use wz_m, only : w3m1
  use mesh, only : if_ortho
  implicit none

  integer, parameter :: lxyz=lx1*ly1*lz1
  real(DP) :: ux(lx1,ly1,lz1,*),uy(lx1,ly1,lz1,*),uz(lx1,ly1,lz1,*),u(lxyz,*)

  real(DP) :: ur(lx1,ly1,lz1),us(lx1,ly1,lz1),ut(lx1,ly1,lz1)

  integer :: e, n, nel

  n = nx1-1
  do e=1,nel
    if (if3d) then
      call local_grad3(ur,us,ut,u,n,e,dxm1,dxtm1)
      if (if_ortho) then
        ux(:,:,:,e) = w3m1*(ur*rxm1(:,:,:,e))
        uy(:,:,:,e) = w3m1*(us*sym1(:,:,:,e))
        uz(:,:,:,e) = w3m1*(ut*tzm1(:,:,:,e))
      else
        ux(:,:,:,e) = w3m1*(ur*rxm1(:,:,:,e) + us*sxm1(:,:,:,e) + ut*txm1(:,:,:,e))
        uy(:,:,:,e) = w3m1*(ur*rym1(:,:,:,e) + us*sym1(:,:,:,e) + ut*tym1(:,:,:,e))
        uz(:,:,:,e) = w3m1*(ur*rzm1(:,:,:,e) + us*szm1(:,:,:,e) + ut*tzm1(:,:,:,e))
      endif
    else
#if 0
      if (ifaxis) then
        call setaxdy(ifrzer(e))  ! reset dytm1
        call setaxw1(ifrzer(e))  ! reset w3m1
      endif

      call local_grad2(ur,us,u,n,e,dxm1,dytm1)

      do i=1,lxyz
        ux(i,e) =w3m1(i,1,1)*(ur(i)*rxm1(i,1,1,e) &
        + us(i)*sxm1(i,1,1,e) )
        uy(i,e) =w3m1(i,1,1)*(ur(i)*rym1(i,1,1,e) &
        + us(i)*sym1(i,1,1,e) )
      enddo
#endif
    endif
  enddo

  return
end subroutine wgradm1

!-----------------------------------------------------------------------
subroutine wlaplacian(out,a,diff,ifld)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lelt, nx1, ny1, nz1
  use input, only : iftmsh
  use tstep, only : ifield, nelfld, imesh
  implicit none

  real(DP) :: out(lx1,ly1,lz1,*),a(*),diff(*)
  real(DP), allocatable :: wrk(:,:,:,:)
  real(DP), allocatable :: h2(:,:,:,:)
  integer :: ifld

  integer :: ntot, ifield_

  ntot = nx1*ny1*nz1*nelfld(ifld)
  if ( .NOT. iftmsh(ifld)) imesh = 1
  if (     iftmsh(ifld)) imesh = 2


  allocate(wrk(nx1,ny1,nz1,nelfld(ifld)), h2(lx1,ly1,lz1,lelt))
  h2 = 0._dp

  ifield_ = ifield
  ifield = ifld

  call bcneusc(out,1)
  call axhelm(wrk,a,diff,h2,imesh,1)
  out(:,:,:,1:nelfld(ifld)) = out(:,:,:,1:nelfld(ifld)) - wrk
   
  ifield = ifield_

  return
end subroutine wlaplacian

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