c-----------------------------------------------------------------------
      subroutine uservp (ix,iy,iz,ieg)
      include 'SIZE'
      include 'NEKUSE'          ! UDIFF, UTRANS

      UDIFF =0.
      UTRANS=0.

      return
      end
c-----------------------------------------------------------------------
      subroutine userf  (ix,iy,iz,ieg)
      include 'SIZE'
      include 'NEKUSE'          ! FFX, FFY, FFZ

      integer ix,iy,iz,ieg

      FFX = 0.0
      FFY = 0.0
      FFZ = 0.0

      return
      end
c-----------------------------------------------------------------------
      subroutine userq  (ix,iy,iz,ieg)
      include 'SIZE'
      include 'NEKUSE'          ! QVOL

      QVOL   = 0.0

      return
      end
c-----------------------------------------------------------------------
      subroutine userchk

      include 'SIZE'            ! NX1, NY1, NZ1, NELV, NIO
      include 'INPUT'           ! UPARAM
      include 'TSTEP'           ! ISTEP, IOSTEP, TIME, LASTEP
      include 'SOLN'            ! V[XYZ], V[XYZ]P, PRP, JP, VMULT
      include 'MASS'            ! BM1
      include 'ADJOINT'         ! IFADJ

      integer n, nit
      real Ek(2),timel(2), omega(2), domega
      save EK, timel, omega
      real vtmp(lx1*ly1*lz1*lelv,ldim)
      character*132 restartf

c     set direct/adjoint mode and load restart field
      if (ISTEP.eq.0) then
         if (int(UPARAM(1)).eq.1) then
            if (NIO.eq.0) write(*,*) 'Simulation in adjoint mode'
            IFADJ = .TRUE.
c     keep current base flow
            call opcopy(vtmp(1,1),vtmp(1,2),vtmp(1,ndim),VX,VY,VZ)
c     read the field
c           restartf = 'prtlin_chan_adj0.restart'
c           call load_fld(restartf)
c     copy fileds
            call opcopy(VXP,VYP,VZP,VX,VY,VZ)
            n = NX2*NY2*NZ2*NELV
            call copy(PRP,PR,n)
c     put back base flow
            call opcopy(VX,VY,VZ,vtmp(1,1),vtmp(1,2),vtmp(1,ndim))
            call rzero(PR,n)
         else
            if (NIO.eq.0) write(*,*) 'Simulation in direct mode'
            IFADJ = .FALSE.
c     keep current base flow
            call opcopy(vtmp(1,1),vtmp(1,2),vtmp(1,ndim),VX,VY,VZ)
c     read the field
c           restartf = 'prtlin_chan_dir0.restart'
c           call load_fld(restartf)
c     copy fileds
            call opcopy(VXP,VYP,VZP,VX,VY,VZ)
            n = NX2*NY2*NZ2*NELV
            call copy(PRP,PR,n)
c     put back base flow
            call opcopy(VX,VY,VZ,vtmp(1,1),vtmp(1,2),vtmp(1,ndim))
            call rzero(PR,n)
         endif
      endif

c      get energy
      domega = 1.0
      nit = 10
      if (mod(ISTEP,nit).eq.0) then
         n = NX1*NY1*NZ1*NELV
         Ek(2) = Ek(1)
         timel(2) = timel(1)
         omega(2) = omega(1)
         Ek(1) = 0.5*(glsc3(VXP,VXP,BM1,n)+glsc3(VYP,VYP,BM1,n))
         timel(1) = TIME
c     get growthrate
         if (Ek(2).gt.0.0.and.timel(1).gt.timel(2)) then
            omega(1) = 0.5*log(Ek(1)/Ek(2))/(timel(1)-timel(2))
            domega = abs((omega(1)-omega(2))/omega(1))
         endif
c     set logs
         if (NIO.eq.0.and.ISTEP.gt.2*nit)
     $        write(*,*) 'Energy ',Ek(1),omega(1),domega,TIME
      endif

c     converged field
      if (ISTEP.eq.NSTEPS) then
c     write perturbation field
         call outpost2(VXP,VYP,VZP,PRP,TP,0,'prt')
      endif
      
      return
      end
c-----------------------------------------------------------------------
      subroutine userbc (ix,iy,iz,iside,eg)

      include 'SIZE'
      include 'NEKUSE'          ! UX, UY, UZ, TEMP, X, Y

c     velocity
      UX = 0.0
      UY = 0.0
      UZ = 0.0

      return
      end
c-----------------------------------------------------------------------
      subroutine useric (ix,iy,iz,ieg)

      include 'SIZE'
      include 'NEKUSE'          ! UX, UY, UZ, TEMP, Z
      include 'SOLN'            ! JP

      real amp, ran

c     velocity
c     base flow read from a restart file

c     perturbation
c     random distribution
         amp = 1.0         
         ran = 3.e4*(ieg+X*sin(Y)) - 1.5e3*ix*iy + .5e5*ix 
         ran = 1.e3*sin(ran)
         ran = 1.e3*sin(ran)
         ran = cos(ran)
         UX  = ran*amp
         
         ran = 2.3e4*(ieg+X*sin(Y)) + 2.3e3*ix*iy - 2.e5*ix 
         ran = 1.e3*sin(ran)
         ran = 1.e3*sin(ran)
         ran = cos(ran)
         UY  = ran*amp

         UZ = 0.0

      return
      end
c-----------------------------------------------------------------------
c     This routine to modify element vertices
      subroutine usrdat

      include 'SIZE'
      
      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat2
      include 'SIZE'

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat3
      return
      end
c-----------------------------------------------------------------------