c***************************************************************
c     BEGIN RUNTIME INITIALIZATION
c***************************************************************

	if(ny .eq. 0) then
		num_dims = 2
	else
		num_dims = 3
	endif

      if( myid .eq. 0 ) then
       write(6,*) ' '
       write(6,*) ' '
       write(6,*)'problem_params.h executed via include statement: flow_solve.x'
      endif
 
      t_write_selective=-999999*dt    ! not implemented yet
      if(num_dims .eq. 2) then
       Ly=1.0                         ! 2d calculation, per unit length in y [m]
       nyplanes=1                     ! parameter ny is zero for 2d xz calcs
      else
       nyplanes=ny
      endif

c     Create specialized data structure for MPI data swaps
      call create_MPI_data_type(nx,ny,nz,locnx,locnz,twoslice,subslice,myid,comm)
      call MPI_BARRIER(comm,ierr)

c     Prescribe initial conditions and nondimensionalize variables
      istart_floats=(t_floats_on-t_start)/dt    ! step to begin float tracking
      call initialize(u,v,w,s1,s2,ambient_density,u_cplx,v_cplx,w_cplx,s1_cplx,s2_cplx,
     *                nx,ny,nz,num_dims,efactor,
     *                Lx,Ly,Lz,dx,dy,dz,U0,DGRAD,s1_scale,s2_scale,rho_0,f,g,kappa_1,kappa_2,nu,
     *                bc_flag,ifax,ifay,ifaz,trigx,trigy,trigz,wn,
     *                wnx,wny,wnz,amp,p,T_diss,diss_flag,N,NM1,M,MM1,MM2,step_flag,
     *                Rot,Ri,Re,Pr_1,Pr_2,istart,iend,i_all,i_selective,i_stat,
     *                dt,t_start,t_end,t_write_selective,
     *                t_stat,xt,work(1,1,1,1),work_cplx(1,1,1,8),myid,
     *                numprocs,locnx,locnz,comm,twoslice,subslice,work(1,1,1,8),
     *                netcdf_file,nyplanes,force_cplx,force_flag,
     *                work(1,1,1,2),scalars)


c     Check accuracy of spectral interpolation:
c     (a) reproduce initial u(z) profile at x=y=0
c     i=1  
c     j=1
c     if( myid .eq. 0 ) open(1,file='output/interp_check')
c     do k=1,nz+1

c      dimensionless positions i.e. on gridpoints
c      x0=(i-1.)*dx/Lz
c      y0=(j-1.)*dy/Lz
c      z0=(k-1.)*dz/Lz

c      interpolate to find u on grid points
c      fcs_flag='c'
c      call ds_interp(x0,y0,z0,u_cplx(1,1,1,N),wnx,wny,wnz,locnx,ny,nz,
c    *                fcs_flag,interp_val,myid,comm)
c      if( myid .eq. 0 ) write(1,*) x0*Lz,y0*Lz,z0*Lz,interp_val*U0

c     enddo

      if( myid .eq. 0 ) write(0,*) 'runtime initialization completed'
      if( myid .eq. 0 ) write(6,*) 'runtime initialization completed'
      if( myid .eq. 0 ) close(1)

c***************************************************************
c     END RUNTIME INITIALIZATION
c***************************************************************