c-----------------------------------------------------------------------
c
c     A short example of Rayleigh-Benard convection for large Ra.
c
c     Rayleigh number is set to 10^k, k=5,6,7, and 8 at times 
c
c         t = 0, 0.1, 0.5, and 1.0
c
c     These breakpoints were determined ad-hoc, by simply running a
c     few cases, but appear to provide an effective path to Ra=10^8.
c
c     Moreover, forcing is turned on exponentially as:
c
c         ffy = RaPr*(1-exp(-t/ts)), with ts = 0.05
c
c     The trick with buoyancy-driven flows is that the nonlinearity
c     of the forcing at early times puts a more severe timestep
c     restriction than does the CFL constraint.   Thus, to keep the
c     buoyancy in check, it's necessary to slowly ramp up to the target
c     target value of Ra.  (The same approach is used in experiments.)
c
c
c     For this case, using a positive value of dt is a good idea so
c     that Nek5000 reduces dt automatically to meet the maximum
c     allowed Courant value.
c
c-----------------------------------------------------------------------
      subroutine uservp (ix,iy,iz,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      udiff  = 0
      utrans = 0

      return
      end
c-----------------------------------------------------------------------
      subroutine userf  (ix,iy,iz,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      common /rayleigh_r/ rapr,ta2pr


      stime  = 0.05
      argt   = -time/stime
      scalet = 1-exp(argt)


      buoy = temp*rapr*scalet

      if (if3d) then
         ffx  =   uy*Ta2Pr
         ffy  = - ux*Ta2Pr
         ffz  = buoy
      elseif (ifaxis) then
         ffx  = -buoy
         ffy  =  0.
      else
         ffx  = 0.
         ffy  = buoy
      endif
c     write(6,*) ffy,temp,rapr,'ray',ieg

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

      qvol   = 0.0
      source = 0.0
      return
      end
c-----------------------------------------------------------------------
      subroutine userbc (ix,iy,iz,iside,ieg)
      include 'SIZE'
      include 'TSTEP'
      include 'INPUT'
      include 'NEKUSE'
      common /rayleigh_r/ rapr,ta2pr

      ux=0.
      uy=0.
      uz=0.

      temp=0.  !     Temp = 0 on top, 1 on bottom

      if (if3d) then
         temp = 1-z
      elseif (ifaxis) then  !      domain is on interval x in [-1,0]
         temp = 1.+x
      else                  ! 2D:  domain is on interval y in [0,1]
         temp = 1.-y
      endif


      return
      end
c-----------------------------------------------------------------------
      subroutine useric (ix,iy,iz,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
      integer idum
      save    idum
      data    idum /99/

      ran = 2.e7*(ieg+x*sin(y)) + 1.e6*ix*iy + 1.e7*ix 
      ran = 1.e9*sin(ran)
      ran = 1.e9*sin(ran)
      ran = cos(ran)
      ran = ran1(idum)
      amp = .005

      temp = 1-y + ran*amp*(1-y)*y*x  ! 2D 

      ux=0.0
      uy=0.0
      uz=0.0

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat
      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat3
      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat2
      include 'SIZE'
      include 'TOTAL'

      common /rayleigh_c/ Ra,Ra1,Ra2,Ra3,Prnum,Ta2,Ek1,Ek2,Ek3,ck

      return
      end
c-----------------------------------------------------------------------
      subroutine userchk
      include 'SIZE'
      include 'TOTAL'

c      common /rayleigh_r/ rapr,ta2pr
c      common /rayleigh_c/ Ra,Ra1,Ra2,Ra3,Prnum,Ta2,Ek1,Ek2,Ek3,ck

c      real Ek0,Ek,t12,t23
c      save Ek0,Ek,t12,t23

c      ifxyo = .true.  ! For VisIt
c      if (istep.gt.iostep) ifxyo = .false.

c      ra = 1.e7
c      if (time.gt.0.1) ra = 1.e6
c      if (time.gt.0.5) ra = 1.e7
c      if (time.gt.1.0) ra = 1.e8

c      prnum   = param(2)     ! Pr=1
c      rapr    = ra*prnum
c      ta2pr   = ta2*prnum

c      n    = nx1*ny1*nz1*nelv
c      Ek0  = Ek
c      Ek   = .5*(glsc3(vx,vx,bm1,n)+ glsc3(vy,vy,bm1,n))/volvm1

c     sigma = 1.e-4
c     de    = abs(Ek-Ek0)/dt

c     if (nid.eq.0.and.(mod(istep,10).eq.0.or.istep.lt.100))
c    $   write(6,6) istep,time,Ek,de,ra
c   6 format(i7,1p4e13.5,' ekde')
      
      character*80 filename(9999)

      ntot = nx1*ny1*nz1*nelv

      if (nid.eq.0) then
        open(unit=199,file='file.list',form='formatted',status='old')
        read(199,*) nfiles
        read(199,'(A80)') (filename(i),i=1,nfiles)
        close(199)
      endif

      call bcast(nfiles,isize)
      call bcast(filename,nfiles*80)

      do i = 1,nfiles
        call load_fld(filename(i))
        call hpts()
      enddo

      call exitt
      return
      end
c-----------------------------------------------------------------------
c
c automatically added by makenek
      subroutine usrsetvert(glo_num,nel,nx,ny,nz) ! to modify glo_num
      integer*8 glo_num(1)
      return
      end