C This is a LES turbulent jet flow example using a recyclcing method
C to get a fully developed turbulent inflow boundary condition.
C The mesh was created out of two circular geometries which were
C extruded into the third dimension followed by a merge operation.
C Close to the outflow boundary the flow will be accelerated to avoid  
C solver instabilities due to backflow.
C Instead of a LES SGS model a simple filtering approach is used. 
c
C-----------------------------------------------------------------------
C
C  USER SPECIFIED ROUTINES:
C
C     - boundary conditions
C     - initial conditions
C     - variable properties
C     - forcing function for fluid (f)
C     - forcing function for passive scalar (q)
C     - general purpose routine for checking errors etc.
C
C-----------------------------------------------------------------------
      subroutine uservp (ix,iy,iz,ieg)

C Set user variable properties

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      return
      end

C=======================================================================

      subroutine userf  (ix,iy,iz,ieg)

C Set user forcing function for the momentum

      include 'SIZE'
      include 'TSTEP'
      include 'NEKUSE'

      FFY = 0.0
      FFX = 0.0
      FFZ = 0.0

      return
      end

C=======================================================================

      subroutine userq  (ix,iy,iz,ieg)

C Set user forcing function for scalars

      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      qvol = 0.0

      return
      end

C=======================================================================

      subroutine userchk

      include 'SIZE'
      include 'TOTAL'

      common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal

      common /myoutflow/ d(lx1,ly1,lz1,lelt),m1(lx1*ly1*lz1,lelt)

      integer*8 glo_num(lx1*ly1*lz1*lelt)
      integer gs_bc_hndl
      save gs_bc_hndl
   
      real m1,imass,imassT
      
      rq = 2.0
      uin = 0.0      

c      nface = 2*ndim
c      nxyz  = nx1*ny1*nz1
c      n     = nxyz*nelv
      
      ifxyo = .true.
      if (istep.gt.iostep) ifxyo = .false.
      
      call turb_outflow(d,m1,rq,uin)
      
      return
      end

C=======================================================================

      subroutine userbc (ix,iy,iz,iside,ieg)

C  Set user boundary conditions

      include 'SIZE'
      include 'TSTEP'
      include 'PARALLEL'
      include 'NEKUSE'

      common /inflow/   uin(lx1,ly1,lz1,lelt)
     &                 ,vin(lx1,ly1,lz1,lelt)
     &                 ,win(lx1,ly1,lz1,lelt)


      ux   = 0.             
      uy   = 0.
      uz   = 0.
      flux = 0.
      temp = 0.

      if (ifield .eq. 1) then     
        if(iside.eq.5) then
            ux = 0  
            uy = 0
           eps = 0.01
             if(r.le.0.5)then
                call sstep(0.0,1.0,0.49,0.0035,R,uz)
                uz = uz + uz*eps*(rand(0)-0.5)
c sstep is a smooth hyperbolic tangent step function
             endif
        endif
      elseif (ifield.eq.2) then
          temp = 0.0                                 
      endif 

      return
      end

C=======================================================================

      subroutine useric (ix,iy,iz,ieg)

C   Set initial conditions

      include 'SIZE'
      include 'TSTEP'
      include 'NEKUSE'

      integer idum
      save    idum
      data    idum / 0 /

c     peturbation level 
      eps = 0.2
      
      
      if (r.le.0.5) then                   
        if (ifield .eq. 1) then
           ux = 0.0
           uy = 0.0
           uz = 0.0
        elseif(ifield .eq. 2) then
          temp = 0.0
        endif
      else
        if (ifield .eq. 1) then
           ux = 0.0
           uy = 0.0
           uz = 0.0
        elseif (ifield .eq. 2) then
           temp = 0.0
        endif
      endif

      return
      end

C=======================================================================

      subroutine usrdat
      
      include 'SIZE'
      include 'TOTAL'
        
      param(66) = 4
      param(67) = 4

 
      return
      end

C=======================================================================

      subroutine usrdat2

      include 'SIZE'
      include 'TOTAL'

c     diam = 10.0
c     rad = diam/2.
c     radm = -rad
c     call rescale_x(xm1,radm,rad)
c     call rescale_x(ym1,radm,rad)
c     call rescale_x(zm1,0,20.0)

      return
      end

C=======================================================================
      subroutine usrdat3

      return
      end
C=======================================================================
       subroutine sstep(vbs,vas,sp,epsmr,r,uz)

c      Smooth Step Function
c      Create a smooth jump from vas,vbs (if r is increasing)
c      with a inflection point sp and a smoothness factor epsmr
c      (epsmr->0 discontinous jump)
c      r denotes the independent variable

       implicit none
       real vbs,vas,sp,r,epsmr
       real delta,aid,uz

       delta = (vas-vbs)
       aid =  0.5*delta*(tanh((sp-r)/epsmr) + 1.0) + vbs
       uz = aid !max(1e-15,aid)

       return
       end
C======================================================================
      subroutine turb_outflow(d,m1,rq,uin)

c     . Set div U > 0 in elements with 'O  ' bc.
c
c     . rq is nominally the ratio of Qout/Qin and is typically 1.5
c
c     . uin is normally zero, unless your flow is zero everywhere
c
c     . d and m1 are work arrays of size (lx1,ly1,lz1,lelt), assumed
c     persistant
c
c     This routine may or may not work with multiple outlets --- it has
c     not been tested for this case.
c
c
c     TYPICAL USAGE -- ADD THESE LINES TO userchk() in your .usr file:
c                      (uncommented)
c
c     common /myoutflow/ d(lx1,ly1,lz1,lelt),m1(lx1*ly1*lz1,lelt)
c     real m1
c     rq  = 2.
c     uin = 0.
c     call turb_outflow(d,m1,rq,uin)
c

      include 'SIZE'
      include 'TOTAL'

      real d(lx2,ly2,lz2,lelt),m1(lx1*ly1*lz1,lelt)

      parameter (lw = 3*lx1*ly1*lz1)
      common /ctmp1/ w(lw)

      integer icalld,noutf,e,f
      save    icalld,noutf
      data    icalld,noutf /0,0/

      real ddmax,cso
      save ddmax,cso
      logical ifout

      character*3 b

      n     = nx1*ny1*nz1*nelv
      n2    = nx2*ny2*nz2*nelv
      nxyz  = nx1*ny1*nz1
      nxyz2 = nx2*ny2*nz2

      if (icalld.eq.0) then
         icalld = 1

         b = 'O  '
         call cheap_dist(m1,1,b)

         call rzero (d,n2)
         ddmax = 0.
         noutf = 0

         do e=1,nelv
           ifout = .false.
           do f=1,2*ndim
             if (cbc(f,e,1).eq.b) ifout = .true. ! outflow
             if (cbc(f,e,1).eq.b) noutf = noutf+1
           enddo
           if (ifout) then
            if (lx2.lt.lx1) then ! Map distance function to Gauss
             call maph1_to_l2(d(1,1,1,e),nx2,m1(1,e),nx1,if3d,w,lw)
            else
             call copy(d(1,1,1,e),m1(1,e),nxyz)
            endif
            dmax  = vlmax(m1(1,e),nxyz)
            ddmax = max(ddmax,dmax)
            call rzero(m1(1,e),nxyz) ! mask points at outflow
           else
             call rone (m1(1,e),nxyz)
           endif
         enddo

         ddmax = glamax(ddmax,1)

         do e=1,nelv
           ifout = .false.
           do f=1,2*ndim
             if (cbc(f,e,1).eq.b) ifout = .true. ! outflow
           enddo
           if (ifout) then
              do i=1,nxyz2
                d(i,1,1,e) = (ddmax - d(i,1,1,e))/ddmax
              enddo
           endif
         enddo
         noutf = iglsum(noutf,1)
      endif

      if (noutf.eq.0) return

      if (uin.ne.0) then ! Use user-supplied characteristic velocity
         ubar = uin
      else
         ubar = glsc3(vx,bm1,m1,n)   ! Masked average
         vbar = glsc3(vy,bm1,m1,n)
         wbar = glsc3(vz,bm1,m1,n)
         volu = glsc2(bm1,m1,n)
         ubar = abs(ubar)+abs(vbar)
         if (if3d) ubar = abs(ubar)+abs(wbar)
         ubar = ubar/volu
      endif

      cs = 3*(rq-1.)*(ubar/ddmax)
      if (.not.ifsplit) cs = cs*bd(1)/dt

      do i=1,n2
         usrdiv(i,1,1,1) = cs*(d(i,1,1,1)**2)
      enddo

      return
      end
c-----------------------------------------------------------------