module communicators
  integer,save  :: gcomm,nloc,nproc
!  integer, save :: phase,confg
end module communicators

module domain_confg
  use communicators
  integer, save :: l1s,l1e,l1d,l2s,l2e,l2d
  integer, save :: nl1,nl2
  integer, save :: ng1,ng2,ngd1,ngd2
end module domain_confg

subroutine setup_confg(n1,n2)
  use domain_confg
  integer n1,n2
  
  ng1  = n1   ; ng2  = n2
  ngd1 = ng1  ; ngd2 = ng2
  nl1  = ngd1 ; nl2  = ngd2/nproc 
  
  l1s = 1          ; l1d = ngd1      ; l1e = ng1
  l2s = 1+nl2*nloc ; l2d = l2s+nl2-1 ; l2e = l2d
  return
end subroutine setup_confg

module domain_phase
  use communicators
  integer, save :: l1s,l1e,l1d,l2s,l2e,l2d
  integer, save :: nl1,nl2
  integer, save :: ng1,ng2,ngd1,ngd2
end module domain_phase 

subroutine setup_phase(n1,n2)
  use domain_phase
  integer n1,n2
  
  ng1  = 2*n2  ; ng2  = n1/2
  ngd1 = ng1   ; ngd2 = ng2
  nl1  = ngd1  ; nl2  = ngd2/nproc
  
  l1s = 1          ; l1d = l1s+nl1-1; l1e = l1d
  l2s = 1+nl2*nloc ; l2d = l2s+nl2-1
  if (l2d.lt.ng2) then
     l2e = l2d
  else
     l2e = ng2
  end if
  if (l2s.gt.ng2) then
     l2s = 0; l2e = 0
  end if

  return
end subroutine setup_phase
   

module fft_stuff
  integer, save :: len1,len2,xc,np1,np2
  real*8, allocatable, save :: work(:),trig1(:),trig2(:)
  real*8, save :: scale0,scale1,scale2
end module fft_stuff

subroutine setup_fft()
  use domain_confg
  use fft_stuff
  integer isign,isys,nt
  
  len1 = ng1;         len2 = ng2;
  xc   = len1/2;      np1 = xc/nproc; 
  np2  = ngd2/nproc

  allocate(trig1(8*ng1)) ; trig1 = 0.0
  allocate(trig2(16*ng2)) ; trig2 = 0.0       
  
  allocate(work(8*ngd1*ngd2/nproc)) ; work=0.0
  
  nt     = 1
  isys   = 0
  isign  = 0
  scale0 = 1.0
  scale1 = 1.0/float(len1)
  scale2 = 1.0/float(len2)
!  call hgfft(isys,len1,scale0,work,work,trig1,work,isys)
!  call ggfft(isys,len2,scale0,work,work,trig2,work,isys)
!  call scfft1dui(ng1, trig1)
!  call cfftm1di(ng2,  trig2)
   call dzfft1dui(ng1, trig1)
   call zfftm1di(ng2,  trig2)
!   call dzfftf(ng1, trig1)
!   call dzfftmf(ng2, trig2)
!    call dzfft(isign,len1,scale0,work,work,trig1,work,isys)
!    call zzfftm(isys,len2,nt,scale0,work,len2,work,len2,trig2,work,isys)
!    call zzfft(isys,len2,scale0,work,work,trig2,work,isys)
!    print *,'trig1 = ', trig1
!    print *,'trig2 = ', trig2
  
  return
end subroutine setup_fft

module IO_cham
  integer, save, dimension(2):: arrdim,gdim,lstart,lend,gstart,gend 
  integer, save:: nbytes
end module IO_cham

subroutine IO_cham_confg()
  use domain_confg
  use IO_cham
  
  nbytes    = 8
  gdim(1)   = ngd1; gdim(2)   = ngd2
  arrdim(1) = nl1; arrdim(2) = nl2
  lstart    = 1
  lend      = arrdim
  gstart(1) = l1s; gstart(2) = l2s
  gend(1)   = l1d; gend(2)   = l2d
  
  return
end subroutine IO_cham_confg

subroutine IO_cham_phase()
  use domain_phase
  use IO_cham
  
  nbytes    = 8
  gdim(1)   = ngd1; gdim(2)   = ngd2
  arrdim(1) = nl1; arrdim(2) = nl2
  lstart    = 1
  lend      = arrdim
  gstart(1) = l1s; gstart(2) = l2s
  gend(1)   = l1d; gend(2)   = l2d
  
  return
end subroutine IO_cham_phase

module wave_stuff
  integer, allocatable, save :: ind1(:),ind2(:)
end module wave_stuff

subroutine setup_wave()
  use domain_phase
  use wave_stuff

!
! NOTE the order of indices in phase space is Y X
!

  allocate(ind1(l1s:l1d+4)) ! read as kx
  allocate(ind2(l2s:l2d+4)) ! read as ky
  
  ind1 = 0; ind2=0
  
  call wave_numbers(ind1,ind2)
  
  return
end subroutine setup_wave

module grid_stuff

  real*8, save :: safety,dx,dy,pi,twopi

contains

 function get_np_for_xy(x,y) result(ip)

	use steps_stuff
	use domain_confg
	use communicators

	implicit none

	real*8,intent(in) :: x,y
	integer :: ip

	integer :: ierr
	real*8 :: ys,ye

	! ************************
	! for vertical decomp only for now...
	! ************************
	do ip = 0,nproc-1
		ys = (nl2*ip+0)*dy
		ye = nl2*(ip+1)*dy
!		write(*,*)'in the ip loop',ys,ye
		if ((y>=ys).AND.(y<ye)) &
			exit
	end do

!	if (ip /= nloc) then
!		write(*,*) '[',nloc,'] ERROR: ip /= nloc: x,y = ',x,' ',y,nit
!		call MPI_FINALIZE(ierr)
!		stop
!!	else
!!		write(*,*) 'ip OK', ip, nloc
!	end if

	! test
!	ip = mod(ip+1,nproc)

  end function get_np_for_xy

end module grid_stuff

module physics_stuff
  real*8, save :: Pe, Pc, Re, Ri, Rc
  real*8, save :: Xmx, Ymx
  real*8, save :: amp, omega
  real*8, save :: tiny, crt
  real*8, save :: F_0, k0, delk, C_t,tperiod
  real*8, save :: f_al, al, delal, aldis
  real*8, save :: al_forcing, alpha, beta, kappa, gc
  real*8, save :: masst, mass0
  real*8, save :: sbeta,beta_forcing, u0,fx0,t0
  real*8, save :: max_g, min_g, bg1, bg2, amp_ran,insize,asp
  real*8, save :: la, ME
  integer,save :: bins, nk
  integer falpha
end module physics_stuff

module statistics_stuff
 real*8, allocatable, save :: pdfu(:),specu(:),pdfu_range(:)
 real*8, allocatable, save :: pdfv(:),specv(:),pdfv_range(:)
 real*8, allocatable, save :: pdft(:),spect(:),pdft_range(:)
 real*8, allocatable, save :: pdfp(:),specp(:),pdfp_range(:)
 real*8, allocatable, save :: k_range(:)
 real*8, save :: binsize
 real*8, save :: delta1, delta2, delta
end module statistics_stuff

subroutine setup_statistics
 use statistics_stuff
 use physics_stuff

 allocate(pdfu(1:bins)) ; allocate(pdfv(1:bins)) 
 allocate(pdft(1:bins)) ; allocate(pdfp(1:bins))

 allocate(specu(1:nk))  ; allocate(specv(1:nk))
 allocate(spect(1:nk))  ; allocate(specp(1:nk))

 allocate(pdfu_range(1:bins+1)) ; allocate(pdfv_range(1:bins+1))
 allocate(pdft_range(1:bins+1)) ; allocate(pdfp_range(1:bins+1))

 allocate(k_range(1:nk+1))
  
 print *,'in setup_statistics',bins,nk, binsize
end subroutine setup_statistics
 

module deriv_stuff
  real*8, allocatable, save :: ak1(:),ak2(:)
  real*8, allocatable, save :: sk1(:),sk2(:)
end module deriv_stuff

subroutine setup_deriv()
  use domain_phase
  use wave_stuff
  use deriv_stuff
  use grid_stuff
  use physics_stuff

!
! NOTE the order of indices in phase space is y x
!
  allocate(ak1(l1s:l1e+4)) ! read as ky
  allocate(ak2(l2s:l2e+4)) ! read as kx

  allocate(sk1(l1s:l1e+4)) ! read as ky
  allocate(sk2(l2s:l2e+4)) ! read as kx

  
  ak1 = twopi*float(ind1)/ymx
  ak2 = twopi*float(ind2)/xmx

  
  ak1(l1e+1) = float(ind1(l1e+1))
  ak2(l2e+1) = float(ind2(l2e+1))        

  
  sk1=ak1; sk2=ak2

  if (ind1(l1e+1) .ne. 0) then
     sk1(ind1(l1e+1))=0.0 ; sk1(ind1(l1e+1)+1)=0.0 
  end if
  if (ind2(l2e+1) .ne. 0) sk2(ind2(l2e+1))=0.0
  
  return
end subroutine setup_deriv

module IO_typedef
  integer, save :: phase,confg
end module IO_typedef


subroutine setup_physics(params,n)
  use physics_stuff
  use statistics_stuff
  integer n, tmp
  real*8 pi,params(n)

  pi    = 4.*atan(1.0)
  aldis = params(49)
  Rc    = params(46)
  Ri    = params(50)
  Re    = params(51)
  Pe    = params(52)
  Pc    = params(45)
  falpha = nint(params(53))
  crt   = params(54)

!  Xmx   = params(55) * pi * 2.0d0**(1.5)
!  Xmx   = params(55) * pi * 2.0d0
!  Ymx   = params(56) * pi
  
  Xmx   = params(55) * pi
!  Ymx   = params(56) * pi
  Ymx   = params(56)* pi

  la    = params(57)
  ME    = params(58)

  amp   = params(61)
  omega = params(62)
  tiny  = params(64)

  al_forcing   = params(73)
  beta_forcing = params(74)
  u0           = params(75)
  fx0	       = params(77)
  t0	       = params(78)

  alpha = params(80)
  beta  = params(81)
  kappa = params(82)
  gc    = params(83)
  sbeta = params(84)
  insize = params(85)
  asp   = params(86)

  F_0   = params(24)
  k0    = params(25)
  delk  = params(26)
  f_al  = params(27)
  al    = params(28)
  delal = params(29)
!  C_t   = params(76)
  tperiod=params(76)
  bg1   = params(47)
  bg2   = params(48)

  bins  = params(90)
  nk    = params(91)
  
  binsize= bins
  tmp = params(92); delta1 = float(tmp)
  tmp = params(93); delta2 = float(tmp)
  tmp = params(94); delta  = float(tmp)
  if(nloc.eq.0) print *,'deltas',delta1,delta2,delta
  
  return
end subroutine setup_physics


subroutine setup_grid(params,n)
  use domain_confg
  use grid_stuff
  use physics_stuff
  integer n
  real*8 params(n)
  
  pi = 4.0*atan(1.0)
  twopi = 2.0*pi
  
  dx = Xmx/float(ng1)
  dy = Ymx/float(ng2)
  
  safety = params(63)
  
  return
end subroutine setup_grid

module steps_stuff
  integer, save :: nit,ntotal
end module steps_stuff

subroutine setup_steps(params,n)
  use steps_stuff
  integer n
  real*8 params(n)
  
  nit = 0
  ntotal = nint(params(30))
  
  return
end subroutine setup_steps

module alias_stuff
  real*8, save :: al_x,al_y
  real*8, save :: ak1_mx2, ak2_mx2
  logical, allocatable, save :: al_mask(:,:)
end module alias_stuff

subroutine setup_alias(params,n)

  use domain_phase
  use alias_stuff
  use physics_stuff
  use grid_stuff
  use deriv_stuff
  real*8 params(n)

  integer i,j
  real*8 kx,ky,kxx,kyy,dxx,dyy
  
  al_x = params(65)
  al_y = params(66)

!  print *,'al_x, al_y',al_x,al_y
!  print *,'ng1/4-1, ng2',nl1/4-1, ng2
  ak1_mx2 = (al_y*float(ng1/4 -1) *twopi/ymx)**2 ! y direction
  ak2_mx2 = (al_x*float(ng2   -1) *twopi/xmx)**2 ! x direction
!  print *,'in alias_stuff',ak1_mx2, ak2_mx2

  allocate(al_mask(l1s:l1d,l2s:l2d))  
  al_mask = .true.   ! if al_mask true then element is set to zero
  
!
! Now set the .false. elemsnts, i.e the ones that should be kept.
!

            do j=l2s,l2e
              kx  = ak2(j)
              kxx = kx*kx
              dxx = kxx/ak2_mx2
              do i=l1s,l1e
                ky  = ak1(i)
                kyy = ky*ky
                dyy = kyy/ak1_mx2

                if (dxx+dyy .le. 1.0) then
!                  print *,i,j,'aliasing'
                  al_mask(i,j) = .false.
                end if 

            end do
          end do

  
  return
end subroutine setup_alias

module definitions
  type out_file
     character (len=80) :: name
     character (len=15) :: acc
     character (len=15) :: form
     logical            :: open
     integer            :: unit
     integer            :: recl
     integer            :: nstep
     integer            :: ndump
     integer            :: head
  end type out_file
end module definitions

module output_stuff
  use definitions
  type (out_file), save :: fff(0:9)     
  integer, save :: iword,run_no
end module output_stuff

subroutine setup_output(params,n,f_name)
  use domain_phase
  use output_stuff
  use statistics_stuff

  integer n
  real*8 params(n)
  integer i
  character(len=80) :: f_name,f0 
  
  iword  = nint(params(10))               ! Word length for output
  run_no = nint(params(41))               ! Run number
  do i=0,9
     
     f0 = trim(f_name)//'_'//achar(48+i)//'_'
     call file_name(f0,fff(i)%name,run_no,.false.)
     
     fff(i)%open = .false.
     
  end do

  do i=1,9 
     
     fff(i)%form  = 'unformatted'
     fff(i)%unit  = 20+i
     fff(i)%nstep = nint(params(30+i))
     fff(i)%ndump = 0
     
  end do
  !
  !  Specifics
  !
  
  fff(0)%acc  = 'sequential'
  fff(0)%form = 'formatted'
  
  fff(1)%head = 5
  
  fff(2)%acc  = 'direct'
  fff(2)%recl = 2000
  
  fff(3)%head = 5

  fff(4)%acc  = 'direct'
  fff(4)%recl = binsize+1

  fff(5)%acc  = 'direct'
  fff(5)%recl = binsize+1

  fff(6)%acc  = 'direct'
  fff(6)%recl = binsize+1

  fff(7)%acc  = 'direct'
  fff(7)%recl = binsize+1

  fff(8)%acc  = 'direct'
  fff(8)%recl = binsize+1

  fff(9)%acc  = 'direct'
  fff(9)%recl = binsize+1
  
  return
end subroutine  setup_output    
      
module input_stuff
  use definitions 
  type (out_file),save, dimension(1) :: f_inp
  integer, save :: jword, prev_run_no, dump_no, read_ps
  logical lphase
end module input_stuff

subroutine setup_input(params,n,f_name,mid)
  use domain_phase
  use input_stuff
  integer n,mid
  real*8 params(n)
  integer i0
  character(len=80) :: f_name,f0 
  
  jword  = nint(params(10))           ! word length (same as iword)
  prev_run_no = nint(params(40))      ! Previos run number
  dump_no     = nint(params(42))      ! Dump number within previos run 
  read_ps     = nint(params(44))

  
  if (prev_run_no .eq. 0) then
     
     f_inp(1)%name = 'static'
     
  else
     
     f0 = trim(f_name)//'_1_'
     call file_name(f0,f_inp(1)%name,prev_run_no,.false.)
     f0 = trim(f_inp(1)%name)
     call file_name(f0,f_inp(1)%name,dump_no,.true.)
     f_inp(1)%head = 5     
  end if
  return
end subroutine setup_input   

module time_stuff
  real*8, save ::  time,prev_time, time0
  integer,save ::  gindicator,prev_nit
  real*8, save ::  dt,  dt1,dt2
  real*8, save ::  gdt
  real*8, save ::  cc0,cc1,cc2
end module time_stuff

module diffuse_stuff
  
  real*8, allocatable,save :: &
       &    ext_t0_1(:),ext_t0_2(:), &
       &    ext_t1_1(:),ext_t1_2(:), &      
       &    ext_t2_1(:),ext_t2_2(:), &
       &    ext_t0(:,:),ext_t1(:,:),ext_t2(:,:),&
       &    exs_t0(:,:),exs_t1(:,:),exs_t2(:,:),&
       &    exr_t0(:,:),exr_t1(:,:),exr_t2(:,:)
  
end module diffuse_stuff

subroutine setup_diffuse()
  use domain_phase
  use diffuse_stuff
!----------------------------------------------------- 1st direction 
!  allocate(ext_t0_1(l1s:l1e)) ; ext_t0_1 = 0.0
!  allocate(ext_t1_1(l1s:l1e)) ; ext_t1_1 = 0.0
!  allocate(ext_t2_1(l1s:l1e)) ; ext_t2_1 = 0.0
!----------------------------------------------------- 2nd direction  
!  allocate(ext_t0_2(l2s:l2e)) ; ext_t0_2 = 0.0
!  allocate(ext_t1_2(l2s:l2e)) ; ext_t1_2 = 0.0
!  allocate(ext_t2_2(l2s:l2e)) ; ext_t2_2 = 0.0

  allocate(ext_t0(l1s:l1e,l2s:l2e)); ext_t0 = 0.0
  allocate(ext_t1(l1s:l1e,l2s:l2e)); ext_t1 = 0.0
  allocate(ext_t2(l1s:l1e,l2s:l2e)); ext_t2 = 0.0
  allocate(exs_t0(l1s:l1e,l2s:l2e)); exs_t0 = 0.0
  allocate(exs_t1(l1s:l1e,l2s:l2e)); exs_t1 = 0.0
  allocate(exs_t2(l1s:l1e,l2s:l2e)); exs_t2 = 0.0
  allocate(exr_t0(l1s:l1e,l2s:l2e)); exr_t0 = 0.0
  allocate(exr_t1(l1s:l1e,l2s:l2e)); exr_t1 = 0.0
  allocate(exr_t2(l1s:l1e,l2s:l2e)); exr_t2 = 0.0

  return
end subroutine setup_diffuse

module gband_number
  integer, save :: mygb(1),myab(1),mybb(1),mycb(1)
!  integer, save :: tgb(1), tab(1), tbb(1), tcb(1)
  integer, save :: tgb(1), tab(1), tcb(1)
  real*8, save  :: pdensity0, pdensity
end module gband_number

!MAY NOT BE NEEDED
!module bands_stuff
!  integer,allocatable,save:: acounts(:),bcounts(:),ccounts(:)
!  integer,allocatable,save:: abandx(:),abandy(:)
!  integer,allocatable,save:: bbandx(:),bbandy(:)
!  integer,allocatable,save:: cbandx(:),cbandy(:)
!  integer,allocatable,save::tabandx(:),tabandy(:)
!  integer,allocatable,save::tbbandx(:),tbbandy(:)
!  integer,allocatable,save::tcbandx(:),tcbandy(:)
!end module bands_stuff

module gband_stuff

  real*8,allocatable,save :: gbx(:),  gby(:)
  real*8,allocatable,save :: tgbx(:),tgby(:)
  real*8,allocatable,save :: dd(:),dd1(:),dd2(:)
  real*8,allocatable,save :: tempx1(:),tempy1(:)
  integer,allocatable,save :: tempx2(:),tempy2(:)
  integer,allocatable,save:: indx(:), indy(:)
  integer,allocatable,save:: tindx(:),tindy(:)
  integer,allocatable,save:: dis1(:),dis2(:),dis3(:),dis4(:)
  integer,allocatable,save:: dis(:),gcounts(:)

end module gband_stuff

subroutine setup_gband()
  use communicators
  use gband_number
  use gband_stuff
  use domain_confg

  if(allocated(gcounts)) deallocate(gcounts)
  if(allocated(dis)) deallocate(dis)
  if(allocated(gbx)) deallocate(gbx)
  if(allocated(gby)) deallocate(gby)
  if(allocated(indx))deallocate(indx)
  if(allocated(indy))deallocate(indy)

  allocate(gcounts(nproc)); gcounts = 0
  allocate(dis(nproc));     dis     = 0
  allocate(gbx(mygb(1))); gbx = 0.0
  allocate(gby(mygb(1))); gby = 0.0
  allocate(indx(mygb(1)));indx= 0
  allocate(indy(mygb(1)));indy= 0

  if(allocated(tgbx))  deallocate(tgbx)
  if(allocated(tgby))  deallocate(tgby)
  if(allocated(dd))    deallocate(dd)
  if(allocated(dd1))   deallocate(dd1)
  if(allocated(dd2))   deallocate(dd2)
  if(allocated(tempx1))deallocate(tempx1)
  if(allocated(tempy1))deallocate(tempy1)
  if(allocated(tempx2))deallocate(tempx2)
  if(allocated(tempy2))deallocate(tempy2)
  if(allocated(tindx)) deallocate(tindx)
  if(allocated(tindy)) deallocate(tindy)
  if(allocated(dis1))  deallocate(dis1)
  if(allocated(dis2))  deallocate(dis2)
  if(allocated(dis3))  deallocate(dis3)
  if(allocated(dis4))  deallocate(dis4)

  allocate(tgbx(tgb(1))); tgbx= 0.0
  allocate(tgby(tgb(1))); tgby= 0.0
  allocate(dd(tgb(1)));   dd  = 0.0
  allocate(dd1(tgb(1)));  dd1 = 0.0
  allocate(dd2(tgb(1)));  dd2 = 0.0
  allocate(tempx1(tgb(1))); tempx1=0.0
  allocate(tempy1(tgb(1))); tempy1=0.0
  allocate(tempx2(tgb(1))); tempx2=0
  allocate(tempy2(tgb(1))); tempy2=0
  allocate(tindx(tgb(1)));tindx= 0
  allocate(tindy(tgb(1)));tindy= 0
  allocate(dis1(tgb(1))); dis1 = 0
  allocate(dis2(tgb(1))); dis2 = 0
  allocate(dis3(tgb(1))); dis3 = 0
  allocate(dis4(tgb(1))); dis4 = 0

  return
end subroutine setup_gband

!MAY NOT BE NEEDED
!subroutine setup_bands()
!  use communicators
!  use gband_number
!  use gband_stuff
!  use bands_stuff
!  use domain_confg
!
!  if(allocated(acounts)) deallocate(acounts)
!  if(allocated(bcounts)) deallocate(bcounts)
!  if(allocated(ccounts)) deallocate(ccounts)
!  if(allocated(abandx)) deallocate(abandx)
!  if(allocated(abandy)) deallocate(abandy)
!  if(allocated(bbandx)) deallocate(bbandx)
!  if(allocated(bbandy)) deallocate(bbandy)
!  if(allocated(cbandx)) deallocate(cbandx)
!  if(allocated(cbandy)) deallocate(cbandy)
!  if(allocated(tabandx)) deallocate(tabandx)
!  if(allocated(tabandy)) deallocate(tabandy)
!  if(allocated(tbbandx)) deallocate(tbbandx)
!  if(allocated(tbbandy)) deallocate(tbbandy)
!  if(allocated(tcbandx)) deallocate(tcbandx)
!  if(allocated(tcbandy)) deallocate(tcbandy)
!
!  allocate(acounts(nproc)); acounts = 0
!  allocate(bcounts(nproc)); bcounts = 0
!  allocate(ccounts(nproc)); ccounts = 0
!
!  allocate(abandx(myab(1)));abandx= 0
!  allocate(abandy(myab(1)));abandy= 0
!  allocate(bbandx(mybb(1)));bbandx= 0
!  allocate(bbandy(mybb(1)));bbandy= 0
!  allocate(cbandx(mycb(1)));cbandx= 0
!  allocate(cbandy(mycb(1)));cbandy= 0
!  allocate(tabandx(tab(1)));tabandx= 0
!  allocate(tabandy(tab(1)));tabandy= 0
!  allocate(tbbandx(tbb(1)));tbbandx= 0
!  allocate(tbbandy(tbb(1)));tbbandy= 0
!  allocate(tcbandx(tcb(1)));tcbandx= 0
!  allocate(tcbandy(tcb(1)));tcbandy= 0
!
!  return
!end subroutine setup_bands

module particle_number
  integer, save :: mypp(2),mymp(2),totpp(2),totmp(2)
  integer, save :: num_pp, num_mp
  integer, save :: myep(1),totep(1)
  integer, save :: lpp(2),lnp(2),gpp(2),gnp(2)
  integer, save :: tnumpp,tnumnp
end module particle_number

module p_definitions
  type YParticle
      real*8 :: oldx
      real*8 :: oldy
      real*8 :: newx
      real*8 :: newy
      real*8 :: radius
      real*8 :: a(3)
      real*8 :: l(3)
      real*8 :: strain
!      integer :: oldi
!      integer :: oldj
!      integer :: newi
!      integer :: newj
      integer :: sign
      integer :: ip
      integer :: escape
      integer :: delete
  end type YParticle

	integer,parameter :: N_PASSED_INTS = 3
	integer,parameter :: N_PASSED_REALS = 12

contains

  subroutine copy_p(psource,pdest)

	implicit none

        type (YParticle),intent(in) :: psource
        type (YParticle),intent(out) :: pdest

	pdest%oldx = psource%oldx
	pdest%oldy = psource%oldy
	pdest%newx = psource%newx
	pdest%newy = psource%newy
	pdest%radius = psource%radius
	pdest%l(1:3) = psource%l(1:3)
	pdest%a(1:3) = psource%a(1:3)
	pdest%strain    = psource%strain
!	pdest%oldi = psource%oldi
!	pdest%oldj = psource%oldj
!	pdest%newi = psource%newi
!	pdest%newj = psource%newj
	pdest%sign = psource%sign
	pdest%escape = psource%escape
	pdest%delete = psource%delete
	pdest%ip = psource%ip

  end subroutine copy_p

	subroutine pack_p_reals(p,buf)

		implicit none

		type(YParticle),intent(in) :: p
		real*8,intent(out) :: buf(N_PASSED_REALS)

		buf(1) = p%newx
		buf(2) = p%newy
		buf(3) = p%oldx
	        buf(4) = p%oldy
		buf(5) = p%radius
		buf(6) = p%l(1)
		buf(7) = p%l(2)
		buf(8) = p%l(3)
		buf(9) = p%a(1)
		buf(10)= p%a(2)
		buf(11)= p%a(3)
		buf(12)= p%strain

	end subroutine pack_p_reals

	subroutine unpack_p_reals(buf,p)

		implicit none

		real*8,intent(in) :: buf(N_PASSED_REALS)
		type(YParticle),intent(inout) :: p

		p%newx   = buf(1)
		p%newy   = buf(2)
	        p%oldx   = buf(3)
		p%oldy   = buf(4)
		p%radius = buf(5)
		p%l(1)   = buf(6)
		p%l(2)   = buf(7)
		p%l(3)   = buf(8)
		p%a(1)   = buf(9)
		p%a(2)   = buf(10)
		p%a(3)   = buf(11)
		p%strain = buf(12)

	end subroutine unpack_p_reals

	subroutine pack_p_ints(p,buf)

		type(YParticle),intent(in) :: p
		integer,intent(out) :: buf(N_PASSED_INTS)

		buf(1) = p%sign
		buf(2) = p%escape
		buf(3) = p%delete

	end subroutine pack_p_ints
	
	subroutine unpack_p_ints(buf,p)

		integer,intent(in) :: buf(N_PASSED_INTS)
		type(YParticle),intent(inout) :: p

		 p%sign   = buf(1)
		 p%escape = buf(2)
		 p%delete = buf(3)

	end subroutine unpack_p_ints

end module p_definitions

module particle

  use p_definitions
  type (YParticle), pointer :: lparticle(:)
  integer,parameter :: P_MEMORY_INC = 200
  integer :: np_last = 0
  integer :: np_max = 0
  integer, save:: tnumpp, tnummp
  integer,allocatable,save :: ppindx(:),ppindy(:)
  integer,allocatable,save :: mpindx(:),mpindy(:)    
  integer,allocatable,save :: tppindx(:),tppindy(:)
  integer,allocatable,save :: tmpindx(:),tmpindy(:)
  integer,allocatable,save :: ppcounts(:),mpcounts(:),tpcounts(:)
  real*8,allocatable, save :: ppx(:),ppy(:)
  real*8,allocatable, save :: mpx(:),mpy(:)
  real*8,allocatable, save :: pwk1(:),pwk2(:),pwk3(:),pwk4(:)
  real*8,allocatable, save :: pwk5(:),pwk6(:),pwk7(:),pwk8(:)
  real*8,allocatable, save :: pwk9(:),pwk10(:)
  real*8,allocatable, save :: ewk1(:),ewk2(:),ewk3(:),ewk4(:)
  real*8,allocatable, save :: ewk5(:),ewk6(:),ewk7(:),ewk8(:)
  real*8,allocatable, save :: mwk1x(:),mwk1y(:),mwk2x(:),mwk2y(:)
  real*8,allocatable, save :: rk3u1(:),rk3u2(:),rk3v1(:),rk3v2(:)
  real*8,allocatable, save :: rk3x1(:),rk3x2(:),rk3y1(:),rk3y2(:)
  real*8,allocatable, save :: myescp(:),tescp(:)
  real*8, save:: bmax,bmin

contains

  function get_new_p() result(index)

        implicit none

        integer :: index
        type (YParticle), pointer :: lparticle_new(:)
        integer :: ip

        if (np_last == np_max) then

                ! the user has requested a new particle and there are no more allocated, so
                ! make more...
                np_max = np_max + P_MEMORY_INC
                allocate(lparticle_new(np_max))

                ! copy all particle data from old to new...
                do ip = 1,np_last
                        call copy_p(lparticle(ip),lparticle_new(ip))
                end do

                if (associated(lparticle)) deallocate(lparticle)
                lparticle => lparticle_new

        end if

        np_last = np_last + 1
        index = np_last
        
  end function get_new_p

end module particle

subroutine setup_output_particle(tnumber)

  use communicators
  use particle
  use particle_number
  use gband_stuff

  integer tnumber

  if(allocated(pwk1)) deallocate(pwk1)
  if(allocated(pwk2)) deallocate(pwk2)
  if(allocated(pwk3)) deallocate(pwk3)
  if(allocated(pwk4)) deallocate(pwk4)
  if(allocated(pwk5)) deallocate(pwk5)
  if(allocated(pwk6)) deallocate(pwk6)
  if(allocated(pwk7)) deallocate(pwk7)
  if(allocated(pwk8)) deallocate(pwk8)
  if(allocated(pwk9)) deallocate(pwk9)
  if(allocated(pwk10))deallocate(pwk10)

  allocate(pwk1(mypp(2))); pwk1 = 0.0
  allocate(pwk2(mypp(2))); pwk2 = 0.0
  allocate(pwk3(totpp(2))); pwk3 = 0.0
  allocate(pwk4(totpp(2))); pwk4 = 0.0
  allocate(pwk5(mymp(2))); pwk5 = 0.0
  allocate(pwk6(mymp(2))); pwk6 = 0.0
  allocate(pwk7(totmp(2))); pwk7 = 0.0
  allocate(pwk8(totmp(2))); pwk8 = 0.0
  if(tnumber .gt. 0) then
    allocate(pwk9(tnumber)); pwk9 = 0.0
    allocate(pwk10(tnumber));pwk10= 0.0 
  endif

  return
end subroutine setup_output_particle

subroutine setup_output_gparticle
 
  use communicators
  use particle
  use particle_number
  use gband_stuff

  if(allocated(pwk1)) deallocate(pwk1)
  if(allocated(pwk2)) deallocate(pwk2)
  if(allocated(pwk3)) deallocate(pwk3)
  if(allocated(pwk4)) deallocate(pwk4)
  if(allocated(pwk5)) deallocate(pwk5)
  if(allocated(pwk6)) deallocate(pwk6)
  if(allocated(pwk7)) deallocate(pwk7)
  if(allocated(pwk8)) deallocate(pwk8)

  allocate(pwk1(mypp(2)+mymp(2))); pwk1 = 0.0
  allocate(pwk2(mypp(2)+mymp(2))); pwk2 = 0.0
  allocate(pwk3(totpp(2))); pwk3 = 0.0
  allocate(pwk4(totpp(2))); pwk4 = 0.0
  allocate(pwk5(mypp(2)+mymp(2))); pwk5 = 0.0
  allocate(pwk6(mypp(2)+mymp(2))); pwk6 = 0.0
  allocate(pwk7(totmp(2))); pwk7 = 0.0
  allocate(pwk8(totmp(2))); pwk8 = 0.0

  return
end subroutine setup_output_gparticle

subroutine setup_escape_particle
 
  use communicators
  use particle
  use particle_number
  use gband_stuff

  if(allocated(ewk1)) deallocate(ewk1)
  if(allocated(ewk2)) deallocate(ewk2)
  if(allocated(ewk3)) deallocate(ewk3)
  if(allocated(ewk4)) deallocate(ewk4)
  if(allocated(ewk5)) deallocate(ewk5)
  if(allocated(ewk6)) deallocate(ewk6)
  if(allocated(ewk7)) deallocate(ewk7)
  if(allocated(ewk8)) deallocate(ewk8)

  allocate(ewk1(myep(1))); ewk1 = 0.0
  allocate(ewk2(myep(1))); ewk2 = 0.0
  allocate(ewk3(myep(1))); ewk3 = 0.0
  allocate(ewk4(myep(1))); ewk4 = 0.0
  allocate(ewk5(totep(1))); ewk5 = 0.0
  allocate(ewk6(totep(1))); ewk6 = 0.0
  allocate(ewk7(totep(1))); ewk7 = 0.0
  allocate(ewk8(totep(1))); ewk8 = 0.0


  return
end subroutine setup_escape_particle

subroutine setup_particle
    
  use communicators
  use particle
  use particle_number
  use gband_stuff

  num_pp = 05
  num_mp = 05
  !if(allocated(ppx)) deallocate(ppx)
  !if(allocated(ppy)) deallocate(ppy)
  !if(allocated(mpx)) deallocate(mpx)
  !if(allocated(mpy)) deallocate(mpy)
  !if(allocated(tppy)) deallocate(tppy)
  !if(allocated(tmpx)) deallocate(tmpx)
  !if(allocated(tmpy)) deallocate(tmpy)

  if(allocated(mwk1x)) deallocate(mwk1x)
  if(allocated(mwk1y)) deallocate(mwk1y)
  if(allocated(mwk2x)) deallocate(mwk2x)
  if(allocated(mwk2y)) deallocate(mwk2y)
  if(allocated(myescp))  deallocate(myescp)
  if(allocated(tescp)) deallocate(tescp)


  if(allocated(ppcounts)) deallocate(ppcounts)
  if(allocated(mpcounts)) deallocate(mpcounts)
  if(allocated(tpcounts)) deallocate(tpcounts)
  if(allocated(ppindx)) deallocate(ppindx)
  if(allocated(ppindy)) deallocate(ppindy)
  if(allocated(mpindx)) deallocate(mpindx)
  if(allocated(mpindy)) deallocate(mpindy)
  if(allocated(tppindx))deallocate(tppindx)
  if(allocated(tppindy))deallocate(tppindy)
  if(allocated(tmpindx))deallocate(tmpindx)
  if(allocated(tmpindy))deallocate(tmpindy)

  !allocate(ppx(mypp(2))); ppx = 0.0
  !allocate(ppy(mypp(2))); ppy = 0.0
  !allocate(mpx(mymp(2))); mpx = 0.0
  !allocate(mpy(mymp(2))); mpy = 0.0
  !allocate(tmpx(totmp(2))); tmpx = 0.0
  !allocate(tmpy(totmp(2))); tmpy = 0.0

  allocate(mwk1x(mypp(1))); mwk1x = 0.0
  allocate(mwk1y(mypp(1))); mwk1y = 0.0
  allocate(mwk2x(mymp(1))); mwk2x = 0.0
  allocate(mwk2y(mymp(1))); mwk2y = 0.0
  allocate(tescp(totep(1))); tescp = 0.0

  allocate(ppcounts(nproc));    ppcounts=0
  allocate(mpcounts(nproc));    mpcounts=0
  allocate(tpcounts(nproc));    tpcounts=0
  allocate(ppindx(mypp(1)));    ppindx = 0
  allocate(ppindy(mypp(1)));    ppindy = 0
  allocate(mpindx(mymp(1)));    mpindx = 0
  allocate(mpindy(mymp(1)));    mpindy = 0
  allocate(tppindx(totpp(1))); tppindx = 0
  allocate(tppindy(totpp(1))); tppindy = 0
  allocate(tmpindx(totmp(1))); tmpindx = 0
  allocate(tmpindy(totmp(1))); tmpindy = 0
 
  return

end subroutine setup_particle 

subroutine setup_particle1

  use communicators
  use particle
  use particle_number
  use gband_stuff


  if(allocated(rk3u1)) deallocate(rk3u1)
  if(allocated(rk3u2)) deallocate(rk3u2)
  if(allocated(rk3v1)) deallocate(rk3v1)
  if(allocated(rk3v2)) deallocate(rk3v2)

  if(allocated(rk3x1)) deallocate(rk3x1)
  if(allocated(rk3x2)) deallocate(rk3x2)
  if(allocated(rk3y1)) deallocate(rk3y1)
  if(allocated(rk3y2)) deallocate(rk3y2)


  allocate(rk3u1(mypp(2)+mymp(2))); rk3u1 = 0.
  allocate(rk3u2(mypp(2)+mymp(2))); rk3u2 = 0.
  allocate(rk3v1(mypp(2)+mymp(2))); rk3v1 = 0.
  allocate(rk3v2(mypp(2)+mymp(2))); rk3v2 = 0.

  allocate(rk3x1(mypp(2)+mymp(2))); rk3x1 = 0.
  allocate(rk3y1(mypp(2)+mymp(2))); rk3y1 = 0.
  allocate(rk3x2(mypp(2)+mymp(2))); rk3x2 = 0.
  allocate(rk3y2(mypp(2)+mymp(2))); rk3y2 = 0.

  return

end subroutine setup_particle1
!*************************************************
subroutine file_name(f1,f2,n,dot)
!-----------------------------------------------------------
!  This routine is used to generate numbered filenames.
! 
!  if f1='banana' and n=12 then on exit
!     f2='banana0012'
!-----------------------------------------------------------
  logical dot
!
  character*1 name(4)
  integer nn(4),i,i0,n
  character(len=80) :: f1,f2
!-------------------------
  nn(1) =  n/1000
  nn(2) = (n - nn(1)*1000)/100
  nn(3) = (n - nn(1)*1000 - nn(2)*100)/10
  nn(4) = (n - nn(1)*1000 - nn(2)*100 - nn(3)*10)
!
  do i=1,4
     name(i) = char(48+nn(i))
  enddo
!    
  i0=index(f1,' ')-1 
  if (dot) then
     f2 = f1(1:i0)//'.'//name(1)//name(2)//name(3)//name(4)
  else 
     f2 = f1(1:i0)//name(1)//name(2)//name(3)//name(4)
  end if
  !
  return
end subroutine file_name