c***********************************************************************
      PROGRAM Surfer
c***
      INCLUDE 'undefined.h'
      INCLUDE 'dimension.h'
      INCLUDE 'global.h'
      INCLUDE 'flags.h'
      INCLUDE 'physics.h'
      INCLUDE 'pointers.h'
      INCLUDE 'mpif.h'
      INTEGER t,size,frame,tswap,a
      INTEGER IERR00,DIMPGR(3),LNPROCS,LPPRANK,LCOMM2D
      INTEGER GP_IFETCH
      DOUBLE PRECISION c
      LOGICAL PERIOD(1),RORDER
      PARAMETER (RORDER = .FALSE.)
      INTEGER v12345678912345a
      INTEGER v12345678912345b

C.... Command line stuff

      integer count
      CHARACTER(len=132) arg

C.... Fortran/C interface

      interface
         subroutine gp_parfile(string) bind(C, name="GP_parfile")
            use iso_c_binding, only: c_char
            character(kind=c_char) :: string(*)
         end subroutine gp_parfile
       end interface
C      
      count = command_argument_count()
      if (count.lt.1) then
         CALL get_command_argument(0, arg)
         print *,'Warning: ',TRIM(arg),' called with no arguments'
         print *,'will use standard input as parameter file.'
      else
         CALL get_command_argument(1, arg)
         print *,'Debugging: parfile name is ',TRIM(arg)
         call gp_parfile(TRIM(arg))
      endif


C.... Starting MPI.
      CALL MPI_INIT(IERR00)
      DEBUG_MPI = GP_IFETCH("debug_mpi",0,"DEFAULT")
      CALL CHECK_MPI(IERR00,11)
C.... Preliminaries

      v12345678912345a=1
      v12345678912345b=2
      if(v12345678912345a.eq.v12345678912345b) STOP 'no long variables'



      PERIOD(1) = .FALSE.
C.... Creating the virtual topology.

c      write(*,*) "Creating the virtual topology"
      CALL MPI_COMM_SIZE(MPI_COMM_WORLD,LNPROCS,IERR00)
      CALL CHECK_MPI(IERR00,12)
      NPROCS = LNPROCS
c        write(*,*) "mpi_comm_world,nprocs,ierr00",
c     &     mpi_comm_world,nprocs,ierr00
c        write(*,*) "after comm size IERR00 = ",IERR00
c      
      DIMPGR(1) = LNPROCS
      CALL MPI_CART_CREATE(MPI_COMM_WORLD,1,DIMPGR,PERIOD,RORDER,
     &                     LCOMM2D,IERR00)
      CALL CHECK_MPI(IERR00,13)

      COMM2D = LCOMM2D
c      write(*,*) "dimpgr(1),nprocs,ierr00,comm2D"
c     &     ,dimpgr(1),lnprocs,ierr00,comm2D
c      write(*,*) "Creating the virtual topology: COMM_RANK"
      CALL MPI_COMM_RANK(LCOMM2D,LPPRANK,IERR00)
      CALL CHECK_MPI(IERR00,14)

c      write(*,*)" pprank = ",lpprank
      PPRANK=LPPRANK
      CALL MPI_BARRIER(LCOMM2D,IERR00)
      CALL CHECK_MPI(IERR00,15)

C.... Obtaining information about virtual topology.
      CALL INFOTP()
      CALL MPI_BARRIER(LCOMM2D,IERR00)
      CALL CHECK_MPI(IERR00,16)


c      CALL MPI_FINALIZE(IERR00)
C      stop "debug"
c***
c     SUBROUTINES: initialize, timestep
c***
      call initialize()
      size = nx*ny*nz
      frame = 0

      write (*,*) 'PPRANK',PPRANK,' nx ny nz t0 :',nx,ny,nz,t0

      tswap = 0
      do T = T0,TMAX
         CALL MPI_BARRIER(LCOMM2D,IERR00)
         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,pprank)
              write (*,*) ' --- Writing dump.dat at t = ',t,' ---'
            endif 
         endif
      enddo
c***
      CALL MPI_FINALIZE(IERR00)
      CALL CHECK_MPI(IERR00,17)

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
      DOUBLE PRECISION FINDNINF, scale
      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')
      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 (GP_ifetch('ring',0,'WARNING').eq.1) then
            call prepdistring(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 (*,*) '--- PPRANK',PPRANK,' RESTARTING WITH FILE dump.',
     &    PPRANK,' .dat ---'
        call iniget(u,v,w,zcc(ia(ng)),zu(ia(ng)),nx,ny,nz,t0,pprank)
	t0=t0+1
        write (*,*) '--- PPRANK',PPRANK,' RESTARTING AT T =',t0,
     &' WITH FILE dump.',PPRANK,' .dat ---'
      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,*) 'PPRANK',PPRANK,' 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
