C Linearized Model problem described in Ananias G. Tomboulides et. al
C JCP 146, CP986079, 1998
C
C Validation case for Low-Mach number implementation
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

      implicit none
      
      integer ix,iy,iz,ieg

      INCLUDE 'SIZE'   
      INCLUDE 'TOTAL'  
      INCLUDE 'NEKUSE' 

      UDIFF = 1.
      UTRANS= 1./TEMP

      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE USERF  (IX,IY,IZ,ieg)                

      implicit none

      integer ix,iy,iz,ieg

      INCLUDE 'SIZE'   
      INCLUDE 'TSTEP'  
      INCLUDE 'NEKUSE' 

      FFY = 0.0
      FFX = 0.0
      FFZ = 0.0

      RETURN
      END
C-----------------------------------------------------------------------

      SUBROUTINE USERQ  (IX,IY,IZ,ieg)                

      implicit none

      integer ix,iy,iz,ieg

      INCLUDE 'SIZE'   
      INCLUDE 'TOTAL'  
      INCLUDE 'NEKUSE' 

      real delta,sech,XD
      parameter(delta=0.2)

      XD = Y/delta
      sech = 1.0/cosh(XD)
     
      qvol = sech*sech/delta * (0.5 + tanh(XD)/delta)

      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE USERCHK

      implicit none

      INCLUDE 'SIZE'
      INCLUDE 'TOTAL'

      COMMON /SCRNS/ w1     (LX1,LY1,LZ1,LELT)
     $              ,T_err  (LX1,LY1,LZ1,LELT)
     $              ,VX_err (LX1,LY1,LZ1,LELT)
     $              ,QTL_err (LX1,LY1,LZ1,LELT)
      real w1,T_err,VX_err,QTL_err

      integer i,j,k,iel,ntot

      real glsc3,glmax

      real sumqw,delta,l2_err,exact,max_err,XD,sech
      parameter(delta=0.2)

      ntot = nx1*ny1*nz1*nelv

      if(istep.lt.1) return

      do 100 iel=1,nelv
      do 100 k=1,nz1
      do 100 j=1,ny1
      do 100 i=1,nx1
         XD = XM1(i,j,k,iel)/delta
         sech = 1.0/cosh(XD)
         exact = 0.5*(3.0 + tanh(XD))
         VX_err(i,j,k,iel) = abs(VX (i,j,k,iel  ) - exact) 
         T_err(i,j,k,iel) = abs( T (i,j,k,iel,1) - exact) 
         QTL_err(i,j,k,iel) = abs(QTL(i,j,k,iel) - 
     &                            0.5/delta*(1. - (tanh(XD)*tanh(XD))))
 100  continue 
  
      max_err = glmax(VX_err,ntot)
      l2_err  = glsc3(VX_err,bm1,VX_err,ntot)/VOLVM1
      if (l2_err.gt.0) l2_err = sqrt(l2_err)
      if(nid.eq.0)
     & write(6,1) istep,time,max_err,l2_err
    1 format(i8,1p3e12.4,' ERROR VX: MAX/L2')

      max_err = glmax(T_err,ntot)
      l2_err  = glsc3(T_err,bm1,T_err,ntot)/VOLVM1
      if (l2_err.gt.0) l2_err = sqrt(l2_err)
      if(nid.eq.0)
     & write(6,2) istep,time,max_err,l2_err
    2 format(i8,1p3e12.4,' ERROR T: MAX/L2')

      max_err = glmax(QTL_err,ntot)
      l2_err  = glsc3(QTL_err,bm1,QTL_err,ntot)/VOLVM1
      if (l2_err.gt.0) l2_err = sqrt(l2_err)
      if(nid.eq.0)
     & write(6,3) istep,time,max_err,l2_err
    3 format(i8,1p3e12.4,' ERROR QTL: MAX/L2')


      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE USERBC (IX,IY,IZ,ISIDE,ieg)                

      implicit none

      integer ix,iy,iz,iside,ieg

      INCLUDE 'SIZE'
      INCLUDE 'TSTEP'
      INCLUDE 'NEKUSE'

      REAL xd,delta
      parameter(delta=0.2)

      XD = Y/delta
 
      IF (IFIELD .eq. 1) THEN          ! velocity
     
            UX = 0.5
            UY = 0.0

      ELSEIF (IFIELD .eq. 2) THEN      ! Temperature

          TEMP = 300

      ENDIF

      return
      END   
C-----------------------------------------------------------------------
      SUBROUTINE USERIC (IX,IY,IZ,ieg)

C   Set initial conditions

      implicit none

      integer ix,iy,iz,ieg

      INCLUDE 'SIZE'
      INCLUDE 'TSTEP'
      INCLUDE 'NEKUSE'

      real delta,xd
      parameter(delta=0.2)

      XD = Y/delta

      IF (IFIELD.EQ.1) THEN          ! ---------- Velocity
          UX = 0.5
          UY = 0.0
      ENDIF

      IF (IFIELD.EQ.2)               ! ---------- Temperature
     &    TEMP = 300

      RETURN
      END
C-----------------------------------------------------------------------
      subroutine usrdat

      return
      end
C-----------------------------------------------------------------------
      subroutine usrdat2

      return
      end
C-----------------------------------------------------------------------
      subroutine usrdat3

      return
      end