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)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
C
      udiff =0.
      utrans=0.
      return
      end
c-----------------------------------------------------------------------
      subroutine userf  (ix,iy,iz,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
c
      common /lesmodela/ smaVD(lx1,ly1,lz1,lelt,3,3)             
      common /lesmodelr/ delta(lx1*ly1*lz1*lelt),gamma,sma_constant,
     $                   Aplus,Ra,Ri_c
cc
      ifield = 1
C
      Ra = 5.e6
c
      xd = 2
      td = 100*dt

      ramp = ramp_to_one (x,xd)*ramp_to_one(time,td)
c
      ie  = gllel(ieg)     ! local (to this processor) element number

 
cc      ffx = ramp*smaVD(ix,iy,iz,ie,1,1)
cc      ffy = ramp*smaVD(ix,iy,iz,ie,2,1) - Ra*temp
cc      ffz = ramp*smaVD(ix,iy,iz,ie,3,1)
      ffx = smaVD(ix,iy,iz,ie,1,1)
      ffy = smaVD(ix,iy,iz,ie,2,1) - Ra*temp
      ffz = smaVD(ix,iy,iz,ie,3,1)
c
c      if ((time.eq.(2*dt)).and. (ie.eq.20).and.(ix.eq.1).and.
c     $ (iy.eq.1)) then
c      if (istep.eq.10) then
c         write(6,6) 'stresses   ', smaVD(ix,iy,iz,ie,1,1),
c     $                             smaVD(ix,iy,iz,ie,2,1),
c     $                             smaVD(ix,iy,iz,ie,3,1)
c 6       format(a11, 3f36.30)
c         write(6,7) 'ramp, temp', ramp, temp
c 7       format(a10, 2f36.34)
c      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine userq  (ix,iy,iz,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
c
      real kappa0
c
      common /lesmodelb/ smaVD_T(lx1,ly1,lz1,lelt,1,3)             
c
C
      ifield = 2
c
      qvol   = 0.0
      source = 0.0
c
      sfy = 1.
      if (y.gt.0.8) sfy = 0.
c
      kappa0 = 100.0
c
      xd = 2
      td = 100*dt

      ramp = ramp_to_one (x,xd)*ramp_to_one(time,td)
c
      ie   = gllel(ieg)     ! local (to this processor) element number

cc      qvol = ramp*smaVD_T(ix,iy,iz,ie,1,1)
cc     $     -exp(-x)*kappa0*(temp-sfy)
      qvol = smaVD_T(ix,iy,iz,ie,1,1)
     $       -exp(-x)*kappa0*(temp-sfy)

c      if ((time.eq.(2*dt)).and. (ie.eq.20).and.(ix.eq.1).and.
c     $ (iy.eq.1)) then
c      if (istep.eq.200) then
c         write(6,7) 'stress_T   ', smaVD_T(ix,iy,iz,ie,1,1)
c 7       format(a11, f36.30)
c      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine userf_set
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
c
      integer  icalld,istep_last
      save     icalld,istep_last
      data     icalld,istep_last   / 0,-1 /
c
      common /lesmodela/ smaVD(lx1,ly1,lz1,lelt,3,3)             
      common /lesmodelr/ delta(lx1*ly1*lz1*lelt),gamma,sma_constant,
     $                   Aplus,Ra,Ri_c

c
      if (icalld.eq.0) then
         icalld=1
c
         ntot = nx1*ny1*nz1*nelv
         N2 = nx1/2
c
         dx = xm1(n2+1,1,1,1)-xm1(n2,1,1,1)
         dx = glmin(dx,1)
         dz = zm1(1,1,n2+1,1)-zm1(1,1,n2,1)
         dz = glmin(dz,1)
c
         call filter_radius(delta)
cTO DO!!!
c         one  = 1.
c         arg  = one/3.
c         do i=1,ntot
c            delta(i) = (abs(dx*dz*delta(i)))**arg
c            write(6,9) 'delta  ', delta(i)
c 9          format(a7, f36.34)
c         enddo
         gamma = 6.
         Aplus = 1.
         sma_constant = 0.05
	 Ri_c=0.25
      endif
c
      if (istep.ne.istep_last) then
         istep_last = istep
        call compute_sma    (smaVD)
        call div_sma_t      (smaVD)
      endif
c

      return
      end
c-----------------------------------------------------------------------
      subroutine userq_set
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
c
      real kappa0
c
      integer  istep_last
      save     istep_last
      data     istep_last   / -1 /
c
      common /lesmodelb/ smaVD_T(lx1,ly1,lz1,lelt,1,3)             
c
      if (istep.ne.istep_last) then
         istep_last = istep
        call compute_sma_T    (smaVD_T)
        call div_sma_t_T      (smaVD_T)
      endif
C

      return
      end
c-----------------------------------------------------------------------
      subroutine userchk
      include 'SIZE'
      include 'TOTAL'
c
      common /bcin/ velx(lx1,ly1,lz1,lelt),unorm
      common /lesmodelr/ delta(lx1*ly1*lz1*lelt),gamma,sma_constant,
     $                   Aplus,Ra,Ri_c
c
      if (istep.eq.1) then  ! Modify the vertical diffusivity
         call reset_vert_diff
         call outpost(vx,vy,vz,pr,delta,'   ')
      endif
c
      call userf_set   ! Set up LES arrays for next step
      call userq_set
c
      ntot  = nx1*ny1*nz1*nelv
      unorm = 0.
      if (istep.ge.1) unorm = glamax(vx,ntot)
c
      return
      end
c-----------------------------------------------------------------------
      subroutine userbc (ix,iy,iz,iside,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
c
      common /bcin/ velx(lx1,ly1,lz1,lelt),unorm
      integer icalld
      save    icalld
      data    icalld /0/
c
      if (icalld.eq.0) then
         icalld = 1
         ntot = nx1*ny1*nz1*nelv
         call set_bc_in
         do i=1,ntot
            velx(i,1,1,1) = bc_in(ym1(i,1,1,1))
         enddo
      endif
c
      ie = gllel(ieg)
      ux=0.
      uy=0.0
      uz=0.0
c
      ymax = 1.
c      ymin = .8
      ymin = .6
      eta  = (ymax-y)/(ymax-ymin)
      temp = 0.5*(1.-cos(pi*eta))
      ux   = 0.33*unorm*velx(ix,iy,iz,ie)
c
      return
      end
c-----------------------------------------------------------------------
      subroutine useric (ix,iy,iz,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
      ux=0.0
      uy=0.0
      uz=0.0
      temp=0
      ymax = 1.
c      ymin = .8
      ymin = .6
      eta  = (ymax-y)/(ymax-ymin)
      temp = 0.5*exp(-x**20.)*(1.-cos(pi*eta))
c
      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat
      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat2
c
c     Modify geometry
c
      include 'SIZE'
      include 'TOTAL'
c
      ntot = nx1*ny1*nz1*nelv
c
      do i=1,ntot
         beta = (10. - xm1(i,1,1,1))/10.
         alph = (1.0 - ym1(i,1,1,1))
c         ym1(i,1,1,1) = ym1(i,1,1,1) + alph*beta*0.8
         ym1(i,1,1,1) = ym1(i,1,1,1) + alph*beta*0.6
      enddo
c
      return
      end
c-----------------------------------------------------------------------
      subroutine set_bc_in
      common /bc_in_dat/ ymax,ymin,cmax,cmin,cscale
c
c
      ymax = 1.0
c      ymin = 0.8
      ymin = 0.6
c
      my = 100
      dy = (ymax-ymin)/my
      de = 1./my
c
      cmin = 0.
      cmax = 0.
      do j=0,my
         y    = ymin + j*dy
         eta = 1.-j*de
         phij = phi_func(eta)
         cmin = min(cmin,phij)
         cmax = max(cmax,phij)
      enddo
c
      cscale = 1. / cmin
c
      return
      end
c-----------------------------------------------------------------------
      function bc_in(y)
c
      common /bc_in_dat/ ymax,ymin,cmax,cmin,cscale
c
c     
      eta   = (ymax-y)/(ymax-ymin)
      bc_in = cscale*phi_func(eta)
c
      return
      end
c-----------------------------------------------------------------------
      function phi_func(eta)
c
c     Compute P4 - P2
c
      common /phi_stuff/ PNN(0:4)
c
      call legendre_poly(pnn,eta,4)
      phi_func = pnn(4)-pnn(2)
c
      return
      end
c-----------------------------------------------------------------------
      subroutine reset_vert_diff
c
      INCLUDE 'SIZE'
      INCLUDE 'GEOM'
      INCLUDE 'INPUT'
      INCLUDE 'MASS'
      INCLUDE 'TSTEP'
      INCLUDE 'WZ'
c
      COMMON /FASTMD/ IFDFRM(LELT), IFFAST(LELT), IFH2, IFSOLV
      LOGICAL         IFDFRM      , IFFAST      , IFH2, IFSOLV
c
      common /scrles/ rvm1(lx1,ly1,lz1,lelt)  ! Make vertical diff. small
     $             , svm1(lx1,ly1,lz1,lelt)
     $             , tvm1(lx1,ly1,lz1,lelt)
     $             , wj  (lx1,ly1,lz1,lelt)
C
      NXYZ1 = NX1*NY1*NZ1
      NTOT1 = NXYZ1*NELT
c
      vert_diff_ratio = 1.e-4  ! Ratio of vertical to horizontal diffusivity
c
      vdrs = sqrt(vert_diff_ratio)
c
      call cmult2(rvm1,rym1,vdrs,ntot1)  ! Vertical direction is Y
      call cmult2(svm1,sym1,vdrs,ntot1)
      call cmult2(tvm1,tym1,vdrs,ntot1)
c
c     Re-compute geometric factors for integrated del-squared operator.
c
      CALL INVERS2 (WJ,JACM1,NTOT1)
      IF (NDIM.EQ.2) THEN
         CALL VDOT2 (G1M1,RXM1,RYM1,RXM1,RYM1,NTOT1)
         CALL VDOT2 (G2M1,SXM1,svm1,SXM1,svm1,NTOT1)
         CALL VDOT2 (G4M1,RXM1,rvm1,SXM1,svm1,NTOT1)
         CALL COL2  (G1M1,WJ,NTOT1)
         CALL COL2  (G2M1,WJ,NTOT1)
         CALL COL2  (G4M1,WJ,NTOT1)
         CALL RZERO (G3M1,NTOT1)
         CALL RZERO (G5M1,NTOT1)
         CALL RZERO (G6M1,NTOT1)
      ELSE
         CALL VDOT3 (G1M1,RXM1,rvm1,RZM1,RXM1,rvm1,RZM1,NTOT1)
         CALL VDOT3 (G2M1,SXM1,svm1,SZM1,SXM1,svm1,SZM1,NTOT1)
         CALL VDOT3 (G3M1,TXM1,tvm1,TZM1,TXM1,tvm1,TZM1,NTOT1)
         CALL VDOT3 (G4M1,RXM1,rvm1,RZM1,SXM1,svm1,SZM1,NTOT1)
         CALL VDOT3 (G5M1,RXM1,rvm1,RZM1,TXM1,tvm1,TZM1,NTOT1)
         CALL VDOT3 (G6M1,SXM1,svm1,SZM1,TXM1,tvm1,TZM1,NTOT1)
         CALL COL2  (G1M1,WJ,NTOT1)
         CALL COL2  (G2M1,WJ,NTOT1)
         CALL COL2  (G3M1,WJ,NTOT1)
         CALL COL2  (G4M1,WJ,NTOT1)
         CALL COL2  (G5M1,WJ,NTOT1)
         CALL COL2  (G6M1,WJ,NTOT1)
      ENDIF
C
C     Multiply the geometric factors GiM1,i=1,5 with the
C     weights on mesh M1.
C
      DO 580 IEL=1,NELT
         IF (IFAXIS) CALL SETAXW1 ( IFRZER(IEL) )
            CALL COL2 (G1M1(1,1,1,IEL),W3M1,NXYZ1)
            CALL COL2 (G2M1(1,1,1,IEL),W3M1,NXYZ1)
            CALL COL2 (G4M1(1,1,1,IEL),W3M1,NXYZ1)
         IF (NDIM.EQ.3) THEN
            CALL COL2 (G3M1(1,1,1,IEL),W3M1,NXYZ1)
            CALL COL2 (G5M1(1,1,1,IEL),W3M1,NXYZ1)
            CALL COL2 (G6M1(1,1,1,IEL),W3M1,NXYZ1)
         ENDIF
  580 CONTINUE
c
      do ie=1,nelt
         ifdfrm(ie) = .true.
      enddo
C
      RETURN
      END
c-----------------------------------------------------------------------
      subroutine filter_radius(delta)
c     this is the new one 

c     returns the radius of the filter at a meshpoint with 
c     coordinates (x,y,z)
c      
      INCLUDE 'SIZE'  
      INCLUDE 'TOTAL' 
c
      real delta(lx1,ly1,lz1,lelv)
c
c     Use the ARCTAN function to have delta go to zero smoothly near the slope
c
c
      C_pi         = 3.1415926535897932385d0
c      C_ mild      = 20.0
c      C_aggressive = 100.0
      C_ mild      = 500.0
      delta_target = 0.01

      do ie=1,nelv
      do iz=1,nx1
      do iy=1,ny1
      do ix=1,nx1
        y_tilde = 0.06*(10.0-xm1(ix,iy,iz,ie))
        y_0     = 1.0
        r_0     = 1.0 - y_tilde
        y       = ym1(ix,iy,iz,ie)
	delta(ix,iy,iz,ie) 
     $     = 2.0 * delta_target
     $   * (atan(C_mild * (r_0*r_0 - (y-y_0)*(y-y_0))) / C_pi)

c        write(6,9) 'delta', xm1(ix,iy,iz,ie), y, delta(ix,iy,iz,ie)
c 9      format(a5, 3f18.15)
      enddo
      enddo
      enddo
      enddo
c
c
      RETURN
      END
c-----------------------------------------------------------------------
      subroutine grad_1(dux,duy,duz,u,ie)
      INCLUDE 'SIZE'   
      INCLUDE 'TOTAL'  
c
      real dux(lx1*ly1*lz1),duy(lx1*ly1*lz1),duz(lx1*ly1*lz1)
     $    ,  u(lx1,ly1,lz1)
c
c
c     Compute the elemental gradient of u
c
      common /scrmg/ dudr(lx1,ly1,lz1)
     $             , duds(lx1,ly1,lz1)
     $             , dudt(lx1,ly1,lz1)
C
      COMMON /FASTMD/ IFDFRM(LELT), IFFAST(LELT), IFH2, IFSOLV
      LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
C
C     Store the inverse jacobian to speed this operation up
C
      COMMON /DUDXYJ/ JACMI(LX1,LY1,LZ1,LELT)
      REAL JACMI
C
      REAL  DRST(LX1,LY1,LZ1)
C
      INTEGER ICALLD
      SAVE    ICALLD
      DATA    ICALLD /0/
C
      nxy1  = nx1*ny1
      nyz1  = ny1*nz1
      nxyz1 = nx1*ny1*nz1
      ntot  = nxyz1*nelt
C
      IF (ICALLD.EQ.0) THEN
         NTOTT = NXYZ1*NELT
         CALL INVERS2(JACMI,jacm1,NTOT)
         ICALLD=1
      ENDIF
C
      IF (IFAXIS) CALL SETAXDY (IFRZER(ie) )
C
      if (if3d) then
         call mxm (dxm1 ,nx1,u,nx1,dudr,nyz1)
         do iz=1,nz1
            call mxm (u(1,1,iz),nx1,dytm1,ny1,duds(1,1,iz),ny1)
         enddo
         call mxm (u,nxy1,dztm1,nz1,dudt,nz1)
c
         do i=1,nxyz1
            dux(i) = (dudr(i,1,1)*rxm1(i,1,1,ie)
     $             +  duds(i,1,1)*sxm1(i,1,1,ie)
     $             +  dudt(i,1,1)*txm1(i,1,1,ie))*jacmi(i,1,1,ie)
            duy(i) = (dudr(i,1,1)*rym1(i,1,1,ie)
     $             +  duds(i,1,1)*sym1(i,1,1,ie)
     $             +  dudt(i,1,1)*tym1(i,1,1,ie))*jacmi(i,1,1,ie)
            duz(i) = (dudr(i,1,1)*rzm1(i,1,1,ie)
     $             +  duds(i,1,1)*szm1(i,1,1,ie)
     $             +  dudt(i,1,1)*tzm1(i,1,1,ie))*jacmi(i,1,1,ie)
         enddo
      else
         call mxm (dxm1 ,nx1,u,nx1,dudr,nyz1)
         call mxm (u,nx1,dytm1,ny1,duds,ny1)
c
         do i=1,nxyz1
            dux(i) = (dudr(i,1,1)*rxm1(i,1,1,ie)
     $             +  duds(i,1,1)*sxm1(i,1,1,ie))*jacmi(i,1,1,ie)
            duy(i) = (dudr(i,1,1)*rym1(i,1,1,ie)
     $             +  duds(i,1,1)*sym1(i,1,1,ie))*jacmi(i,1,1,ie)
         enddo
         call rzero(duz,nx1*ny1)
      endif
c
      return
      end

c-----------------------------------------------------------------------

      subroutine compute_sma (smaVD)

c     Computes the Smagorinsky term WITH VAN DRIEST DAMPING.
c     We will add this term to the RHS of the equations, in USERF.
c     Note that we do not split this term into two parts anymore.
c
c     d/dj (2.0 * C_s * delta^2 |\grad_sym u| * 0.5 * (du_j/dx_i + du_i/dx_j) )
c                 
c     where \grad_sym u   = 1/2 (du_i/dx_j + du_j/dx_i) and
c           |\grad_sym u| = (2.0 * \grad_sym u * \grad_sym u)^{1/2} 
c     (Hughes' definition)
c
c
      INCLUDE 'SIZE'   
      INCLUDE 'TOTAL' 

      real smaVD (lx1,ly1,lz1,lelt,3,3)
c
      common /scrles/ w(lx1,ly1,lz1,3,3),w_T(lx1,ly1,lz1,1,3)
      common /lesmodelr/ delta(lx1,ly1,lz1,lelt),gamma,sma_constant,
     $                   Aplus,Ra,Ri_c
      real norm_grad_sym, s, yplus, VD_const
      real Ri, factor_Ri, factor
c

      nxyz = nx1*ny1*nz1
      ntot = nx1*ny1*nz1*nelv
c
      do j=1,3
       do i=1,3
         call rzero(smaVD(1,1,1,1,i,j),ntot)
       enddo
      enddo
c

      do ie=1,nelt
         call grad_1(w(1,1,1,1,1),w(1,1,1,1,2),w(1,1,1,1,3),
     $               vx(1,1,1,ie),ie)
         call grad_1(w(1,1,1,2,1),w(1,1,1,2,2),w(1,1,1,2,3),
     $               vy(1,1,1,ie),ie)
         if (if3d) 
     $   call grad_1(w(1,1,1,3,1),w(1,1,1,3,2),w(1,1,1,3,3),
     $               vz(1,1,1,ie),ie)
c        for the temperature
	 call grad_1(w_T(1,1,1,1,1),w_T(1,1,1,1,2),w_T(1,1,1,1,3),
     $               T(1,1,1,ie,1),ie)

         do kz=1,nz1	 
         do ky=1,ny1
         do kx=1,nx1
          norm_grad_sym = 0.0
          do i=1,ndim
           do j=1,ndim
            norm_grad_sym = norm_grad_sym + 
     :                  0.25*(w(kx,ky,kz,i,j)+w(kx,ky,kz,j,i))*
     :                       (w(kx,ky,kz,i,j)+w(kx,ky,kz,j,i))
           enddo
          enddo
cc        this is the definition in Hughes
          norm_grad_sym = 2.0*norm_grad_sym 
          if (norm_grad_sym.gt.0) norm_grad_sym = sqrt(norm_grad_sym)
cc        Compute the Richardson number          
c          if ((w(kx,ky,kz,1,2)**2 + w(kx,ky,kz,3,2)**2).gt.1.e-6) then
          if ((w(kx,ky,kz,1,2)**2 + w(kx,ky,kz,3,2)**2).gt.0.) then
c             Ri = Ra*dabs(w_T(kx,ky,kz,1,2)) 
             Ri = -Ra*w_T(kx,ky,kz,1,2) 
     $          / (w(kx,ky,kz,1,2)**2 + w(kx,ky,kz,3,2)**2)
          else
             Ri = Ri_c
          endif

          if(Ri.ge.Ri_c) then 
             factor = 0.0
          elseif(Ri.le.0.) then 
             factor = 1.0
          else
             factor = sqrt(1-Ri/Ri_c)            
          endif

          do j=1,ndim
           do i=1,ndim
            smaVD(kx,ky,kz,ie,i,j) = norm_grad_sym*0.5
c     $                             * factor 
     $                             * (w(kx,ky,kz,j,i)+w(kx,ky,kz,i,j))
           enddo
          enddo
         enddo
         enddo
         enddo
      enddo
c
c      write(6,7) 'Ri_max is   ', Ri_max
c 7    format(a12, g)
c
      s = sma_constant
cc    this is the definition in Hughes
      s = 2.0*s


      do j=1,3
      do i=1,3
      do ie=1,nelt
      do ky=1,ny1
      do kx=1,nx1
       y_tilde = 0.06*(10.0-xm1(kx,ky,1,ie))
c       VD_const = 1.0 - exp(ym1(kx,ky,1,ie)-y_tilde)
c       VD_const = 1.0 - exp(-(ym1(kx,ky,1,ie)-y_tilde)/Aplus) 
       do kz=1,nz1	 
          smaVD(kx,ky,kz,ie,i,j) = smaVD(kx,ky,kz,ie,i,j)
     :             *delta(kx,ky,kz,ie)*delta(kx,ky,kz,ie)
c     :             *s*VD_const*VD_const
     :             *s
       enddo
      enddo
      enddo
      enddo
      enddo
      enddo


      return
      end


c-----------------------------------------------------------------------

      subroutine compute_sma_T (smaVD_T)


c     Same as above, but for the temperature.
c     
c
      INCLUDE 'SIZE'   
      INCLUDE 'TOTAL' 

      real smaVD_T (lx1,ly1,lz1,lelt,1,3)
c
      common /scrles/ w(lx1,ly1,lz1,3,3),w_T(lx1,ly1,lz1,1,3)
      common /lesmodelr/ delta(lx1,ly1,lz1,lelv),gamma,sma_constant,
     $                   Aplus,Ra,Ri_c
      real norm_grad_sym, s, yplus, VD_const
      real Ri, factor_Ri, factor
c

      nxyz = nx1*ny1*nz1
      ntot = nx1*ny1*nz1*nelv
c
      do j=1,3
       call rzero(smaVD_T(1,1,1,1,1,j),ntot)
      enddo
c

      do ie=1,nelt
         call grad_1(w(1,1,1,1,1),w(1,1,1,1,2),w(1,1,1,1,3),
     $               vx(1,1,1,ie),ie)
         call grad_1(w(1,1,1,2,1),w(1,1,1,2,2),w(1,1,1,2,3),
     $               vy(1,1,1,ie),ie)
         if (if3d) 
     $   call grad_1(w(1,1,1,3,1),w(1,1,1,3,2),w(1,1,1,3,3),
     $               vz(1,1,1,ie),ie)
c        for the temperature
	 call grad_1(w_T(1,1,1,1,1),w_T(1,1,1,1,2),w_T(1,1,1,1,3),
     $               T(1,1,1,ie,1),ie)

         do kz=1,nz1	 
         do ky=1,ny1
         do kx=1,nx1
          norm_grad_sym = 0.0
          do i=1,ndim
           do j=1,ndim
            norm_grad_sym = norm_grad_sym + 
     :                  0.25*(w(kx,ky,kz,i,j)+w(kx,ky,kz,j,i))*
     :                       (w(kx,ky,kz,i,j)+w(kx,ky,kz,j,i))
           enddo
          enddo
cc        this is the definition in Hughes
          norm_grad_sym = 2.0*norm_grad_sym 
          if (norm_grad_sym.gt.0) norm_grad_sym = sqrt(norm_grad_sym)
cc        Compute the Richardson number          
c          if ((w(kx,ky,kz,1,2)**2 + w(kx,ky,kz,3,2)**2).gt.1.e-6) then
          if ((w(kx,ky,kz,1,2)**2 + w(kx,ky,kz,3,2)**2).gt.0.) then
c             Ri = Ra*dabs(w_T(kx,ky,kz,1,2)) 
             Ri = -Ra*w_T(kx,ky,kz,1,2) 
     $          / (w(kx,ky,kz,1,2)**2 + w(kx,ky,kz,3,2)**2)
          else
             Ri = Ri_c
          endif
          if(Ri.ge.Ri_c) then 
             factor = 0.0
          elseif(Ri.le.0.) then 
             factor = 1.0
          else
             factor = sqrt(1-Ri/Ri_c)            
          endif
          do j=1,ndim
            smaVD_T(kx,ky,kz,ie,1,j) = norm_grad_sym
c     $                               * factor
     $                               * w_T(kx,ky,kz,1,j)
          enddo
         enddo
         enddo
         enddo
      enddo
c

      s = sma_constant
cc    this is the definition in Hughes
      s = 2.0*s

      do j=1,3
      do ie=1,nelt
      do ky=1,ny1
      do kx=1,nx1
       y_tilde = 0.06*(10.0-xm1(kx,ky,1,ie))
c       VD_const = 1.0 - exp(ym1(kx,ky,1,ie)-y_tilde)
c       VD_const = 1.0 - exp(-(ym1(kx,ky,1,ie)-y_tilde)/Aplus) 
       do kz=1,nz1	 
          smaVD_T(kx,ky,kz,ie,1,j) = smaVD_T(kx,ky,kz,ie,1,j)
     :                          *delta(kx,ky,kz,ie)*delta(kx,ky,kz,ie)
c     :                          *s*VD_const*VD_const
     :                          *s
       enddo
      enddo
      enddo
      enddo
      enddo


c
      return
      end


c-----------------------------------------------------------------------
      subroutine div_sma_t (sma_rhs)
c
c                               T
c     This routine computes -div (sma_rhs), with dssum and multiplication
c         -1
c     by B , so that the result is in the form of "ffx" (i.e., ready for
c
c     dssum and B to be applied).
c
c
c
c     Store result in sma_rhs_i1
c
      INCLUDE 'SIZE'   
      INCLUDE 'TOTAL'  
c
      COMMON /DUDXYJ/ JACMI(LX1,LY1,LZ1,LELT)
      REAL JACMI
c
      real sma_rhs(lx1*ly1*lz1*lelt,3,3)
c
      common /scrles/ du(lx1,ly1,lz1,lelt,3)
     $             , dv(lx1,ly1,lz1,lelt,3)
c
c
      nxy1 = nx1*ny1
      nyz1 = ny1*nz1
      ntot = nx1*ny1*nz1*nelv
c
      do i=1,ndim
         call rzero(du(1,1,1,1,i),ntot)
      enddo
c
      do i=1,ndim
c
       if (if3d) then
         do k=1,ntot
            dv(k,1,1,1,1) = ( rxm1(k,1,1,1)*sma_rhs(k,i,1)
     $                    + rym1(k,1,1,1)*sma_rhs(k,i,2)
     $                    + rzm1(k,1,1,1)*sma_rhs(k,i,3) )
     $                    * jacm1(k,1,1,1)*bm1(k,1,1,1)
            dv(k,1,1,1,2) = ( sxm1(k,1,1,1)*sma_rhs(k,i,1)
     $                    + sym1(k,1,1,1)*sma_rhs(k,i,2)
     $                    + szm1(k,1,1,1)*sma_rhs(k,i,3) )
     $                    * jacm1(k,1,1,1)*bm1(k,1,1,1)
            dv(k,1,1,1,3) = ( txm1(k,1,1,1)*sma_rhs(k,i,1)
     $                    + tym1(k,1,1,1)*sma_rhs(k,i,2)
     $                    + tzm1(k,1,1,1)*sma_rhs(k,i,3) )
     $                    * jacm1(k,1,1,1)*bm1(k,1,1,1)
         enddo
         do ie=1,nelv
            call mxm (dxtm1,nx1,dv(1,1,1,ie,1),nx1,du(1,1,1,ie,1),nyz1)
            do iz=1,nz1
             call mxm (dv(1,1,iz,ie,2),nx1,dym1,ny1,du(1,1,iz,ie,2),ny1)
            enddo
            call mxm (dv(1,1,1,ie,3),nxy1,dzm1,nz1,du(1,1,1,ie,3),nz1)
         enddo
         do l=1,ntot
           sma_rhs(l,i,1) = -(du(l,1,1,1,1)+du(l,1,1,1,2)+du(l,1,1,1,3))
         enddo
       else
         do k=1,ntot
            dv(k,1,1,1,1) = ( rxm1(k,1,1,1)*sma_rhs(k,i,1)
     $                    +   rym1(k,1,1,1)*sma_rhs(k,i,2) )
     $                    *   jacmi(k,1,1,1)*bm1(k,1,1,1)
            dv(k,1,1,1,2) = ( sxm1(k,1,1,1)*sma_rhs(k,i,1)
     $                    +   sym1(k,1,1,1)*sma_rhs(k,i,2) )
     $                    *   jacmi(k,1,1,1)*bm1(k,1,1,1)
         enddo
         do ie=1,nelv
            call mxm (dxtm1,nx1,dv(1,1,1,ie,1),nx1,du(1,1,1,ie,1),ny1)
            call mxm (dv(1,1,1,ie,2),nx1,dym1 ,ny1,du(1,1,1,ie,2),ny1)
         enddo
         do k=1,ntot
            sma_rhs(k,i,1) = -(du(k,1,1,1,1)+du(k,1,1,1,2))
         enddo
       endif

      enddo
c                -1
c     dssum and B    :
c
      ifield = 1
      do i = 1,ndim
         call dssum(sma_rhs(1,i,1),nx1,ny1,nz1)
         call col2 (sma_rhs(1,i,1),binvm1,ntot)
      enddo

      return
      end
c
c-----------------------------------------------------------------------
      subroutine div_sma_t_T (sma_rhs)
c
c     div (sma_rhs)
c
      INCLUDE 'SIZE'   
      INCLUDE 'TOTAL'  
c
      COMMON /DUDXYJ/ JACMI(LX1,LY1,LZ1,LELT)
      REAL JACMI
c
      real sma_rhs(lx1*ly1*lz1*lelt,3)
c
      common /scrles/ du(lx1,ly1,lz1,lelt,3)
     $              , dv(lx1,ly1,lz1,lelt,3)
c
c
      nxy1 = nx1*ny1
      nyz1 = ny1*nz1
      ntot = nx1*ny1*nz1*nelv
c
      do i=1,ndim
         call rzero(du(1,1,1,1,i),ntot)
      enddo
c
      if (if3d) then
         do k=1,ntot
            dv(k,1,1,1,1) = ( rxm1(k,1,1,1)*sma_rhs(k,1)
     $                    +   rym1(k,1,1,1)*sma_rhs(k,2)
     $                    +   rzm1(k,1,1,1)*sma_rhs(k,3) )
     $                    *   jacm1(k,1,1,1)*bm1(k,1,1,1)
            dv(k,1,1,1,2) = ( sxm1(k,1,1,1)*sma_rhs(k,1)
     $                    +   sym1(k,1,1,1)*sma_rhs(k,2)
     $                    +   szm1(k,1,1,1)*sma_rhs(k,3) )
     $                    *   jacm1(k,1,1,1)*bm1(k,1,1,1)
            dv(k,1,1,1,3) = ( txm1(k,1,1,1)*sma_rhs(k,1)
     $                    +   tym1(k,1,1,1)*sma_rhs(k,2)
     $                    +   tzm1(k,1,1,1)*sma_rhs(k,3) )
     $                    *   jacm1(k,1,1,1)*bm1(k,1,1,1)
         enddo
         do ie=1,nelv
            call mxm (dxtm1,nx1,dv(1,1,1,ie,1),nx1,du(1,1,1,ie,1),nyz1)
            do iz=1,nz1
             call mxm (dv(1,1,iz,ie,2),nx1,dym1,ny1,du(1,1,iz,ie,2),ny1)
            enddo
             call mxm (dv(1,1,1,ie,3),nxy1,dzm1,nz1,du(1,1,1,ie,3),nz1)
         enddo
         do k=1,ntot
            sma_rhs(k,1) = -(du(k,1,1,1,1)+du(k,1,1,1,2)+du(k,1,1,1,3))
         enddo
      else
         do k=1,ntot
            dv(k,1,1,1,1) = ( rxm1(k,1,1,1)*sma_rhs(k,1)
     $                    +   rym1(k,1,1,1)*sma_rhs(k,2) )
     $                    *   jacmi(k,1,1,1)*bm1(k,1,1,1)
c    $                    *   bm1(k,1,1,1)
            dv(k,1,1,1,2) = ( sxm1(k,1,1,1)*sma_rhs(k,1)
     $                    +   sym1(k,1,1,1)*sma_rhs(k,2) )
     $                    *   jacmi(k,1,1,1)*bm1(k,1,1,1)
c    $                    *   bm1(k,1,1,1)
         enddo
         do ie=1,nelv
            call mxm (dxtm1,nx1,dv(1,1,1,ie,1),nx1,du(1,1,1,ie,1),ny1)
            call mxm (dv(1,1,1,ie,2),nx1,dym1 ,ny1,du(1,1,1,ie,2),ny1)
         enddo
         do k=1,ntot
            sma_rhs(k,1) = -(du(k,1,1,1,1)+du(k,1,1,1,2))
         enddo
      endif
c                -1
c     dssum and B    :
c
      ifield = 2
      call dssum(sma_rhs,nx1,ny1,nz1)
      call col2(sma_rhs,tmask,ntot)
c      call outpost(sma_rhs,sma_rhs,vz,pr,sma_rhs,'   ')
      call col2 (sma_rhs,bintm1,ntot)
c      call outpost(sma_rhs,sma_rhs,vz,pr,sma_rhs,'   ')
c      write(6,*) 'this is binv',ntot
c      call exitt
c
      return
      end
c
c-----------------------------------------------------------------------
      function ramp_to_zero (x,xd)
c
c     For x < 0 , function yields 1
c     For x > xd, function yields 0
c     In between, a sine function is used.
c
      ramp_to_zero = 0.
c
      if (x.lt.0) then
         ramp_to_zero = 1.0
      elseif (x.lt.xd) then
         arg = 0.5*(3.14159)*x/xd
         ramp_to_zero = 1.-sin(arg)
      endif
c
      return
      end
c-----------------------------------------------------------------------
      function ramp_to_one (x,xd)
c     
c     For x < 0 , function yields 0
c     For x > xd, function yields 1
c     In between, a sine function is used.
c
      ramp_to_zero = 0.
c    
      if (x.lt.0) then
         ramp_to_zero = 1.0
      elseif (x.lt.xd) then
         arg = 0.5*(3.14159)*x/xd
         ramp_to_zero = 1.-sin(arg)
      endif
      ramp_to_one = 1-ramp_to_zero
c     
      return
      end
c-----------------------------------------------------------------------