c***********************************************************************
      PROGRAM Surfer
c***
      INCLUDE 'undefined.h'
      INCLUDE 'dimension.h'
      INCLUDE 'global.h'
      INCLUDE 'flags.h'
      INCLUDE 'physics.h'
      INCLUDE 'pointers.h'
      INTEGER t,size,frame, tswap
      INTEGER v12345678912345a
      INTEGER v12345678912345b

      v12345678912345a=1
      v12345678912345b=2
      if(v12345678912345a.eq.v12345678912345b) STOP 'no long variables'
c***
c     SUBROUTINES: initialize, timestep
c***
      call initialize()
      size = nx*ny*nz
      frame = 0

      write (*,*) 't0 =', t0

      tswap = 0
      do t = t0,tmax
         call timestep(t,tswap)
         if(tdump.gt.0) then
            if (mod(t,tdump) .eq. 0.and. t .gt. 0) then
              call save(u,v,w,zcc(ia(ng)),zu(ia(ng)),nx,ny,nz,t)
              write (*,*) ' --- Writing dump.dat at t = ',t,' ---'
            endif 
         endif
      enddo
c***
      stop
      end
c***********************************************************************
      SUBROUTINE initialize()
c***
      INCLUDE 'undefined.h'
      INCLUDE 'dimension.h'
      INCLUDE 'global.h'
      INCLUDE 'flags.h'
      INCLUDE 'physics.h'
      INCLUDE 'pointers.h'
      INTEGER n1, n2, n3, restart,GP_IFETCH
      INTEGER nappe, rayleighjet
      DOUBLE PRECISION scale
      DOUBLE PRECISION FINDNINF
      EXTERNAL FINDNINF,GP_IFETCH
c***
c     SUBROUTINES: read_flags, prepjet, writeminmax, read_physics,
c                  initialize_bc, rescale, bc_scalar, bc_vector,
c                  initpressure, printdata, nappe
c***
      call read_flags(n1,n2,n3)
      nx = n1
      ny = n2
      nz = n3
c***
c     check nx, ny, nz against the maximum values in global.h
c***
      write (*,*)nx,ny,nz
      if (nx .gt. NMX) then
         write (*,*) 'STOP: nx > NMX:', nx, NMX
         STOP
      else 
         write (*,*) 'nx:',nx,' NMX:', NMX
      endif 
      if (ny .gt. NMY) then
         write (*,*) 'STOP: ny > NMY:', ny, NMY
         STOP
      else 
         write (*,*) 'ny:',ny,' NMY:', NMY
      endif 
      if (nz .gt. NMZ) then
         write (*,*) 'STOP: nz > NMZ:', nz, NMZ
         STOP
      else 
         write (*,*) 'nz:',nz,' NMZ:', NMZ
      endif 
      write (*,*)
c***
c     check max dimension of MEMLEN and initialize indices
c***
      call initpressure()
c***
c     read IC and make a few checks
c***
      restart=0
      restart=GP_IFETCH('restart',0,'SEVERE')
      nappe=GP_IFETCH('nappe',0,'WARNING')
      rayleighjet=GP_IFETCH('rayleighjet',0,'WARNING')
      write (*,*) '--- RESTART = ', restart
      if(restart.eq.0) then
        !  call prepjet(u,v,w,zcc(ia(ng)),zu(ia(ng)),S11,nx,ny,nz)
         if(nappe.eq.1)  then
            call prepnappe(u,v,w,zcc(ia(ng)),zu(ia(ng)),S11,nx,ny,nz)
         else if (rayleighjet.eq.1) then
            call preprayjet(u,v,w,zcc(ia(ng)),zu(ia(ng)),S11,nx,ny,nz)
         else
            call prepdistdrop(u,v,w,zcc(ia(ng)),zu(ia(ng)),S11,nx,ny,nz)
         endif
         t0 = 0
      else
        write (*,*) '--- RESTARTING WITH FILE dump.dat ---'
        call iniget(u,v,w,zcc(ia(ng)),zu(ia(ng)),nx,ny,nz,t0)
      endif
      write (*,*)
      call writeminmax('cc',zcc(ia(ng)),nx,ny,nz)
      call writeminmax('u',u,nx,ny,nz)
      call writeminmax('v',v,nx,ny,nz)
      call writeminmax('w',w,nx,ny,nz)
      call writeminmax('p',zu(ia(ng)),nx,ny,nz)
      write (*,*)
c***
      call read_physics()
c***  
c     if t0 = 0 make velocities adimensional and rescale them according
c     to : v --> v1 = v*timescale/lengthsize --> v2 = v1*tau/h;
c     else print restart message
c***  
      if (t0 .eq. 0) then
        if (restart.eq.0) then
          scale = (1.0d0/vscale) * (tau/h)
          call rescale (scale,u,nx*ny*nz) 
          call rescale (scale,v,nx*ny*nz) 
          call rescale (scale,w,nx*ny*nz) 
        endif
        write (6,*) 'Surfer: initial velocity max  = ',
     %        FINDNINF(u,v,w,nx,ny,nz)
        call bc_c(zcc(ia(ng)),nx,ny,nz)
        call bc_vector(u,v,w,nx,ny,nz)
      endif 
      call printdata()
c***
      return 
      end
c***********************************************************************
      SUBROUTINE rescale (scale,a,N)
c***  
      DOUBLE PRECISION scale
      INTEGER N, i
      DOUBLE PRECISION a(N)
c***  
      do i=1,N
         a(i) = a(i) *scale
      enddo 
c***  
      return
      end
