module wei_work
  use sep
  use wei_process
  use wei_type
  use wei_lagrange
  implicit none
  integer,  private, pointer, save :: igood_len(:),ihxbeg(:,:),ihybeg(:,:),ihxend(:,:),ihyend(:,:), do_freq_block(:)
  real, private, save,pointer :: zdepth(:,:)

  !---------------

contains
  !----------------------------------------------------------------
  
  subroutine distribute_init(allow_var)
   character(len=128) :: d_type
   integer :: nz,i,iw,npts,n1,nw
   logical, optional :: allow_var
   REAL  :: loc
   logical :: reduce_h,var
   integer, external :: getch

   var=.true.
   if(present(allow_var)) var=allow_var

   
  if(pro%restart) then
     call restart_status()
  else
    allocate(do_freq_block(size(pro%iw_thread)))
    do i=1,size(pro%iw_thread)
      do_freq_block(i)=i
    end do
  end if 




   d_type="                                                         "
   call from_param("distribute",d_type,"reg")
   if(.not. var) d_type="reg"
   if(d_type(1:3)=="reg") then
     pro%distribute_pattern=0;
     call distribute_reg_nz()
     call distribute_reg_sample()
   else if(d_type(1:3) == "vsf" ) then
     pro%distribute_pattern=2;
     call var_freq_slow()
   else if(d_type(1:2) == "vs" ) then
     pro%distribute_pattern=1;
     call distribute_reg_nz()
     call distribute_var_slow()
   else 
     write(0,*) "distribute pattern",trim(d_type)
     call seperr("unimplemented distribution pattern ")
   end if
  if(pro%verb) write(pro%io_unit,*) "distribution type=",trim(d_type(1:10))," mode=",pro%distribute_pattern




  if(pro%distribute_pattern>-1) then
     
     call from_param("nz_mem",nz,0)
     call from_param("nlagrange_pts",npts,5)
     allocate(pro%filled(pro%az%n))
     

   !in core z buffer
   if(nz >0) then
     allocate(pro%v_i_buf(pro%amx%n,pro%amy%n,pro%ahx%n,pro%ahy%n,nz))
     allocate(pro%mem_buf(nz))
     pro%nmem_buf=nz
   else 
     pro%nmem_buf=0
     allocate(pro%mem_buf(1))
   end if


   !on disk z buffer 
   if(npts-nz >0) then
     pro%ndisk_buf=npts-nz
     pro%use_disk_buf=.true.
     allocate(pro%disk_buf(npts-nz))
   else
     pro%NDISK_BUf=0
     pro%use_disk_buf=.false.
     allocate(pro%disk_buf(1))
   end if

    pro%nw_local=pro%nw_local*pro%nws
      
  do iw=1,size(pro%z_samp)
     nw=size(pro%z_samp(iw)%dz)
     allocate(pro%z_samp(iw)%iloc(npts,pro%az%n))
     allocate(pro%z_samp(iw)%frac(npts,pro%az%n))
     allocate(pro%z_samp(iw)%zloc(nw))
     allocate(pro%z_samp(iw)%status(nw))
     pro%z_samp(iw)%zloc(1)=pro%az%o
     do i=2,nw
       pro%z_samp(iw)%zloc(i)=pro%z_samp(iw)%zloc(i-1)+pro%z_samp(iw)%dz(i-1)
     end do
     call lagrange_init(pro%z_samp(iw)%iloc,pro%z_samp(iw)%frac,pro%z_samp(iw)%zloc)
   end do
end if

  call init_h_reduce()
  call save_pattern()

  end subroutine distribute_init

  

   subroutine restart_status()
    integer :: n1,nredo,i,redo(1000),nrd,ierr
    integer, allocatable :: j(:),tempi(:)
    integer :: npart,ndo
    logical :: self_restart
    integer,external :: sreed

!    call from_param("self_restart",self_restart,.true.)
!    if(self_restart) call create_restart_status()

    !read in the threads not started
    call from_aux("restart_stat","ndo",ndo)
    call from_aux("restart_stat","npart",npart)
    if(npart+ndo==0) call erexit("asked to restart but specified 0 redo_proc")

#ifdef SEP_MPI
if(pro%mpi_ithread==0) then
   ierr=sreed("restart_stat",j,4*npart+ndo*2)
end if
  call MPI_Bcast(j,size(j),MPI_INTEGER, 0, MPI_COMM_WORLD,ierr)
#else
   ierr=sreed("restart_stat",j,4*npart+ndo*2)
#endif

    !read in the threads unfinished
    allocate(pro%write_freq(size(pro%iw_thread))); pro%write_freq=0
    allocate(j(npart+ndo*2))
    deallocate(pro%iw_thread)
    allocate(pro%iw_thread(npart+ndo))
    allocate(do_freq_block(npart+ndo))
    do i=1,ndo
       pro%iw_thread(i)=j(i) 
       do_freq_block(i)=j(i) 
    end do
    do i=1,npart
       pro%iw_thread(i+ndo)=j(i*2-1) 
       do_freq_block(i+ndo)=j(i*2-1) 
       pro%write_freq(j(i*2-1))=j(i*2)
    end do
    deallocate(j)
  end subroutine

  subroutine save_pattern()
#ifdef SEP_MPI
if(pro%mpi_ithread==0) then
#endif
    call to_history("n1",size(pro%iw_thread),"distrib_pattern")
    call srite("distrib_pattern",real(pro%iw_thread),size(pro%iw_thread)*4)
#ifdef SEP_MPI
end if
#endif
  end subroutine


subroutine init_h_reduce()
logical :: reduce_h
integer ::i,nmax,ifreq
logical :: rite_sample
integer, external :: getch
real :: hslope,h0,hmax,z0
rite_sample=.false.

#ifdef SEP_MPI
if(pro%mpi_ithread==0) then
#endif
if(1==getch("sample_h","d",i)) then
  rite_sample=.true.
end if
#ifdef SEP_MPI
end if
#endif

 allocate(igood_len(230))
igood_len=(/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,&
     18, 20, 21, 22, 24, 26, 28, 30, 33, 35, 36, 39, 40, 42, 44,&
     45, 48, 52, 55, 56, 60, 63, 65, 66, 70, 72, 77, 78, 80, 84,&
     88, 90, 91, 99, 104, 105, 110, 112, 117, 120, 126, 130, 132,&
    140, 143, 144, 154, 156, 165, 168, 176, 180, 182, 195, 198, 208,&
    210, 220, 231, 234, 240, 252, 260, 264, 273, 280, 286, 308, 312,&
    315, 330, 336, 360, 364, 385, 390, 396, 420, 429, 440, 455, 462,&
    468, 495, 504, 520, 528, 546, 560, 572, 585, 616, 624, 630, 660,&
    693, 715, 720, 728, 770, 780, 792, 819, 840, 858, 880, 910, 924,&
    936, 990, 1001, 1008, 1040, 1092, 1144, 1155, 1170, 1232, 1260,&
   1287, 1320, 1365, 1386, 1430, 1456, 1540, 1560, 1584, 1638, 1680,&
   1716, 1820, 1848, 1872, 1980, 2002, 2145, 2184, 2288, 2310, 2340,&
   2520, 2574, 2640, 2730, 2772, 2860, 3003, 3080, 3120, 3276, 3432,&
   3465, 3640, 3696, 3960, 4004, 4095, 4290, 4368, 4620, 4680, 5005,&
   5040, 5148, 5460, 5544, 5720, 6006, 6160, 6435, 6552, 6864, 6930,&
   7280, 7920, 8008, 8190, 8580, 9009, 9240, 9360, 10010, 10296,&
  10920, 11088, 11440, 12012, 12870, 13104, 13860, 15015, 16016,&
  16380, 17160, 18018, 18480, 20020, 20592, 21840, 24024, 25740,&
  27720, 30030, 32760, 34320, 36036, 40040, 45045, 48048, 51480,&
  55440, 60060, 65520, 72072, 80080, 90090, 102960, 120120, 144144,&
 180180, 240240, 360360, 720720/)

call from_param("reduce_h",reduce_h,.false.)
call from_param("reduce_h_h0",h0,max(&
  max(abs(pro%ahx%o),abs(pro%ahx%o+pro%ahx%d*(pro%ahx%n-1))),&
  max(abs(pro%ahy%o),abs(pro%ahy%o+pro%ahy%d*(pro%ahy%n-1)))))
call from_param("reduce_h_z0",z0,pro%az%o+pro%az%d*(pro%az%n-1))
call from_param("reduce_h_slope",hslope,0.)
nmax=1;
do i=1,size(pro%iw_first)	
 nmax= max(size(pro%z_samp(i)%dz),nmax)
end do

allocate(pro%ihyb_beg(nmax))
allocate(pro%ihxb_beg(nmax))
allocate(pro%ihyb_end(nmax))
allocate(pro%ihxb_end(nmax))
allocate(ihxbeg(nmax,size(pro%iw_first)))
allocate(ihxend(nmax,size(pro%iw_first)))
allocate(ihybeg(nmax,size(pro%iw_first)))
allocate(ihyend(nmax,size(pro%iw_first)))
allocate(zdepth(nmax,size(pro%iw_first)))
ihxbeg=1
ihybeg=1
ihxend=pro%ahx%n
ihyend=pro%ahy%n

if(rite_sample) then
  call to_history("n1",nmax,"sample_h")
  call to_history("n2",size(pro%iw_first),"sample_h")
  call to_history("n3",3,"sample_h")
end if
zdepth=0
 do ifreq=1,size(pro%iw_first)
  do i=1,size(pro%z_samp(ifreq)%dz)
    if(i>1) zdepth(i,ifreq)=zdepth(i-1,ifreq)+pro%z_samp(ifreq)%dz(i-1)
if (reduce_h) then
    if(zdepth(i,ifreq) > z0) then
      hmax=h0+hslope*(zdepth(i,ifreq)-z0)
      if(hmax < 0) then
         write(0,*) "h dist went zero change hslope, z0, and/or h0"
          write(0,*) "CHECK z,i,ifreq,h ",zdepth(i,ifreq),i,hmax,pro%mpi_ithread
          call erexit("")
       end if
      if(pro%ahy%n>1) then
        call set_limits(pro%ahy,ihybeg(i,ifreq),ihyend(i,ifreq),hmax)
      end if
      if(pro%ahx%n>1) then
        call set_limits(pro%ahx,ihxbeg(i,ifreq),ihxend(i,ifreq),hmax)
end if
      end if
    end if
  end do
 end do

   if(rite_sample) call srite("sample_h",real(zdepth),size(zdepth)*4)
   if(rite_sample) call srite("sample_h",real(ihxbeg),size(ihxbeg)*4)
   if(rite_sample) call srite("sample_h",real(ihxend),size(ihxend)*4)




end subroutine

subroutine set_limits(ax,ibeg,iend,max_dist) 
integer :: ibeg,iend,ileft,iright
real :: mymax,max_dist
type(axis) :: ax
integer :: inext,i 


  mymax=max(abs(ax%o),abs(ax%o+ax%d*(ax%n-1)))
  ibeg=1  
  iend=ax%n
  ileft=0
  iright=0
  if(mymax > max_dist) then
    if(-ax%o > max_dist) then
      ileft=(-ax%o-max_dist)/ax%d
    end if  
    if(ax%o > max_dist) then
      ileft=(ax%o-max_dist)/ax%d
    end if  
    if(-(ax%o+ax%d*(ax%n-1)) > max_dist) then
      iright=-((ax%o+ax%d*(ax%n-1))-max_dist)/ax%d
    end if  
    if((ax%o+ax%d*(ax%n-1)) > max_dist) then
      iright=((ax%o+ax%d*(ax%n-1))-max_dist)/ax%d
    end if  
    inext=next_fast(ax%n-ileft-iright)
    if(inext < ax%n) then
      do while(iend-ibeg+1 /= inext) 
        if(ileft > iright) then
          ileft=ileft-1
          ibeg=ibeg+1
        else
          iright=iright-1
          iend=iend-1
        end if
      end do
     end if  
   end if  
end subroutine set_limits

integer function next_fast(nlow)
integer i,nlow
logical :: found
  i=1; found=.false.
  do while(i < size(igood_len) .and. .not. found)
    if(igood_len(i)< nlow) then
     i=i+1
    else
     found=.true.
    end if
  end do
  next_fast=igood_len(i)
end function


  subroutine setup_var_block(iws,izs)
    integer :: iws,izs
    integer :: j,k,l,i,iz
    logical :: found
    if(pro%distribute_pattern==0)  then
       pro%dzb=pro%az%d
       pro%ihxb_beg=ihxbeg(pro%nzs*(izs-1)+1:pro%nzs*izs,iws)
       pro%ihyb_beg=ihybeg(pro%nzs*(izs-1)+1:pro%nzs*izs,iws)
       pro%ihxb_end=ihxend(pro%nzs*(izs-1)+1:pro%nzs*izs,iws)
       pro%ihyb_end=ihyend(pro%nzs*(izs-1)+1:pro%nzs*izs,iws)
    else
      j=(izs-1)*pro%nzs
      k=min(size(pro%z_samp(iws)%dz),pro%nzs*izs)
      l=k-j
      pro%dzb(1:l)=pro%z_samp(iws)%dz(j+1:k)
      do iz=1,l
       i=1; found=.false.
       do while (i< size(zdepth,1) .and. .not. found )
        if(pro%z_samp(iws)%zloc(iz+j) < zdepth(i,iws)) then
          found=.true.
        else
         i=i+1
        end if
       end do
        pro%ihxb_beg=ihxbeg(i,iws)
        pro%ihxb_end=ihxend(i,iws)
        pro%ihyb_beg=ihybeg(i,iws)
        pro%ihyb_end=ihyend(i,iws)
      end do
    end if
  end subroutine setup_var_block

  subroutine  var_samp_iw(iwb)
    integer :: iwb
    pro%filled=.false.
    pro%mem_buf=0
    pro%disk_buf=0
    pro%z_samp(iwb)%status=0
  end subroutine var_samp_iw

  subroutine distribute_reg_nz() 
    integer :: iw,i
    integer :: isend
    
    !For now we are just going to cycle through the groups
    isend=0
    pro%nw_local=0
    
    if(size(do_freq_block) < pro%mpi_nthread) then
       write(0,"(a,i2,a,i2)") "number of threads=",pro%mpi_nthread,&
       " number frequency groups=",size(do_freq_block)
       call erexit("number of frequencies groups must be greater (reduce nwp)");
    else if(mod(size(do_freq_block), pro%mpi_nthread)/=0) then
       write(0,"(a,i2,a,i2)") "number of threads=",pro%mpi_nthread,&
       " number frequency groups=",size(do_freq_block)
       write(0,*) "warning number of frequency blocks not multiple of np"
    end if

    do iw=1,size(do_freq_block)
       pro%iw_thread(iw)=isend
       if(pro%debug) write(0,"(a,i3,a,i3)") "SENDING  thread block",iw,&
       " to  processor ", pro%iw_thread(iw)
       if(isend==pro%mpi_ithread) pro%nw_local=pro%nw_local+1
       isend=isend+1
       if(isend==pro%mpi_nthread) isend=0
    end do
    
    !Record the first frequency in each group
    allocate(pro%iw_first(pro%nw_local))
    isend=0
    do iw=1,size(do_freq_block)
       if(pro%iw_thread(iw)==pro%mpi_ithread) then
          isend=isend+1
          pro%iw_first(isend)=(do_freq_block(iw)-1)*pro%nws+1
       end if
    end do
    allocate(pro%z_samp(pro%nw_local))
end subroutine distribute_reg_nz
    
subroutine distribute_reg_sample()
  integer :: iw
    do iw=1,pro%nw_local
      allocate(pro%z_samp(iw)%dz(pro%az%n))
      allocate(pro%z_samp(iw)%zloc(pro%az%n))
      pro%z_samp(iw)%dz=pro%az%d
    end do
end subroutine distribute_reg_sample

subroutine create_modified_slowness(S,dz,zloc)
type(weslow_comp) :: S
character(len=128) :: tag_in,tag_out
integer :: sample_per
integer :: n(3)
integer :: ipos
real :: d(3),freq,min_z,max_z
real :: o(3),z,zstep
real :: dz(:),zloc(:)
character(len=128) :: label
real, allocatable :: in_slice(:,:,:),out_slice(:,:)
integer npts,inew
logical :: changed,velocity
integer :: iz,iin,nv_in(2)
real :: fb,fe,min_dz,neg,mys,ov(3),dv(3)
integer :: ib,ie,j,i,jb,je,icur,iread,nv(3)
integer, external :: sseek,sreed,srite
 


tag_in=s%tag_local
tag_out=s%tag_var
velocity=s%velocity

 nv=(/s%amx%n,s%amy%n,s%az%n/)
 ov=(/s%amx%o,s%amy%o,s%az%o/)
 dv=(/s%amx%d,s%amy%d,s%az%d/)



!NOW TIME WRITE OUR NEW VELOCITY/SLOWNESS FIELD TO DISK
!
!
! THIS CAN BE TRICKY....
!
!
! When sample size is smaller than slowness sampling we will
! just do linear interpolation
!
! When it is larger than slowness sampling we will take an
! average  of all of the full samples covered by the jump
! plus fractional portions of other layers
!

allocate(in_slice(nv(1),nv(2),2))
allocate(out_slice(nv(1),nv(2)))
if(0/= sseek(tag_out,0,0)) &
  call seperr("trouble with initial seek on output in variable slow")
if(zloc(2)< zloc(1))  then
 neg=-1.
else 
 neg=1.;
end if

!WHEN ADJOINT THIS SHOULD ALSO FlIP THE SlOWNESS FIELD

write(0,*) pro%mpi_ithread,size(zloc)

do i=1,size(zloc)
  if(neg==1) then
    fb=(zloc(i)-ov(3))/dv(3)+1;  ib=fb; fb=fb-ib
    fe=(zloc(i)+dz(i)-ov(3))/dv(3)+1;ie=fe; fe=fe-ie
  else
    fe=(zloc(i)-ov(3))/dv(3)+1;  ie=fe; fe=fe-ie
    fb=(zloc(i)-dz(i)-ov(3))/dv(3)+1;ib=fb; fb=fb-ib
  end if 


  mys=0.
  jb=ib+1;je=ie+1
  if(ib<1) then
    ib=1; jb=1
  end if
  if(ib>=nv(3))then
   ib=nv(3)-1; jb=nv(3)-1
  end if
  if(ie>=nv(3))then
   ie=nv(3)-1; je=nv(3)-1
  end if

  if(ib-1 .ne. sseek_block(trim(tag_in),ib-1,nv(1)*nv(2)*4,0)) &
   call seperr("trouble with  seek 1 in variable slow construction")
  if( size(in_slice)*4 .ne.sreed(trim(tag_in),in_slice,size(in_slice)*4)) then
    write(0,*) "From tag",trim(tag_in),"seeked",ib-1," of ",nv(1)*nv(2)*4
    write(0,*) zloc(i),ov(3),dz(i)
    call seperr("trouble with reed 1 in variable slow construction")
  end if 
  iin=ib

  !add in the first portion
  out_slice=((1.-fb)*in_slice(:,:,1)+fb*in_slice(:,:,2))*(1.-fb)*dv(3)
  mys=mys+dv(3)*(1.-fb)


  !add in middle suff
  do j=ib+1,ie-1
     if(j>iin+1) then
      if(j-1 .ne. sseek_block(trim(tag_in),j-1,nv(1)*nv(2)*4,0))&
   call seperr("trouble with  seek 2 in variable slow construction")
      if(size(In_slice)*4 .ne. sreed(trim(tag_in),in_slice,size(in_slice)*4)) &
    call seperr("trouble with reed 2 in variable slow construction")
     iin=j
     end if
     out_slice=out_slice+dv(3)*in_slice(:,:,1+j-iin)
     mys=mys+dv(3)
  end do

  !add in the end stuff
  if(ie/=iin) then
    if(ie-1 .ne. sseek_block(trim(tag_in),ie-1,nv(1)*nv(2)*4,0))&
      call seperr("trouble with  seek 3 in variable slow construction")
    if(size(in_slice)*4 .ne. sreed(trim(tag_in),in_slice,size(in_slice)*4)) &
      call seperr("trouble with reed 3 in variable slow construction")

iin=ie
  end if
  out_slice=out_slice+((1.-fe)*in_slice(:,:,1)+fe*in_slice(:,:,2))*(fe)*dv(3)
  mys=mys+dv(3)*fe

  out_slice=out_slice/mys

  if(.not. velocity) out_slice=1./out_slice
  if(size(out_slice)*4 .ne. srite(trim(tag_out),out_slice,size(out_slice)*4))&
   call seperr("trouble with write in variable slow construction")
end do
if(0/= sseek(trim(tag_out),0,0)) &
  call seperr("trouble with final seek of output in varaible slow construct")
deallocate(out_slice,in_slice)


end subroutine


subroutine distribute_var_slow()
real :: min_zsamp,max_zsamp,base_freq,mys
integer :: cycles,nsamp,samples,i
real    ::  dz(10000)

call from_param("min_zsamp",min_zsamp,pro%az%d/3.)
call from_param("max_zsamp",max_zsamp,10000.)
call from_param("cycles",cycles,1000000)
call from_param("samp_per",samples,40)



call from_param("base_freq",base_freq,40.)
nsamp=calc_work(base_freq/6.285714,min_zsamp,max_zsamp,cycles,samples,dz)
  mys=0
do i=1,nsamp
  mys=mys+dz(i)
end do
do i=1,size(pro%z_samp)
		allocate(pro%z_samp(i)%dz(nsamp))
    pro%z_samp(i)%dz(1:nsamp)=dz(1:nsamp)
end do

end subroutine distribute_var_slow

integer function calc_work(freq,min_dz,max_dz,max_cycle,sample_per,dz)
integer :: sample_per
logical :: slowness
integer :: n(3)
integer :: ipos,ithread,max_cycle
real :: d(3),freq,min_z,max_z,nax_dz,dv(3),max_Dz
real :: o(3),z,zstep,ov(3)
real :: dz(:)
integer i,npts,inew
real :: sample_hold(size(dz,1))
logical :: changed
integer :: iz
real :: fb,fe,min_dz,myf,cut,count
integer :: ib,ie,j,jb,je,ifreq,ierr,nv(3)
logical :: cutoff

nv=(/pro%slow%s_s%amx%n,pro%slow%s_s%amy%n,pro%slow%s_s%az%n/)
ov=(/pro%slow%s_s%amx%o,pro%slow%s_s%amy%o,pro%slow%s_s%az%o/)
dv=(/pro%slow%s_s%amx%d,pro%slow%s_s%amy%d,pro%slow%s_s%az%d/)
min_z=pro%slow%s_s%az%o 
max_z=pro%slow%s_s%az%o +(pro%slow%s_s%az%n-1)*pro%slow%s_s%az%d


n=nv
d=dv
o=ov


myf=freq

z=min_z;npts=1;sample_hold(npts)=z
cutoff=.false.
count=0
do while(z< max_z .and. .not. cutoff) 
  ipos=min(n(3),max(1,floor((z-o(3))/d(3)+1.5)))
  !GIVEN THE CURRENT VELOCITY DEPTH LOCATION, CHECK TO SEE HOW
  !LARGE IN DEPTH THIS FREQUENCY VALUE AND ARE SAMPLING
  zstep=(1./pro%max_slow(ipos)/myf)/sample_per

  !check to see if there is a lower velocity with our given depth step
  inew=min(n(3),max(1,floor((z+zstep-o(3))/d(3)+1.5)))
  zstep=max(min_dz,(minval(1./pro%max_slow(ipos:inew))/myf)/sample_per)
  zstep=min(max_dz,zstep)
  z=z+zstep
  npts=npts+1;sample_hold(npts)=z
  count=count + freq/maxval(pro%max_slow(ipos:inew))* zstep
  if(count>max_cycle) cutoff=.true.
end do
!    minval(pro%max_slow(ipos:inew)), &
!    zstep/( freq/minval(pro%max_slow(ipos:inew))),&
!    freq/minval(pro%max_slow(ipos:inew))


!CALCULATE DZ SAMPLING
do i=1,npts-1 
  dz(i)=sample_hold(i+1)-sample_hold(i)
end do

calc_work=npts


end function calc_work 

subroutine var_freq_slow()
integer :: iw,i,k
integer :: isend,nc,ierr
real :: min_samp,max_samp,base_freq,io,lowcost
integer :: cycles,nsamp,samples,itry,inode(1),igroup,j,nw,hxk
real    ::  freq,lowest,myrand,mycost,min_zsamp,max_zsamp
real,allocatable    :: mywork(:),dz(:,:),cost_node(:)
integer, allocatable   :: nz(:)
integer, allocatable   :: low_plan(:),cur_plan(:)
real :: mid
logical,allocatable :: not_distrib(:)
logical :: rite_sample
integer, external :: getch
rite_sample=.false.


#ifdef SEP_MPI
    if(pro%mpi_ithread==0) then
#endif
if(1==getch("sample_z","d",i)) rite_sample=.true.
#ifdef SEP_MPI
end if
#endif

call from_param("min_zsamp",min_zsamp,pro%az%d/3.)
call from_param("max_zsamp",max_zsamp,10000.)
call from_param("cycles",cycles,1000000)
call from_param("samp_per",samples,8)
call from_param("io_cost",io,20.)

isend=0
pro%nw_local=0

nw=size(do_freq_block)
allocate(dz(100000,nw))
allocate(nz(nw))
allocate(not_distrib(nw))
allocate(low_plan(nw))
allocate(cur_plan(nw))
allocate(mywork(nw))
allocate(cost_node(pro%mpi_nthread))

do iw=1,nw
  mid=pro%nws*(do_freq_block(iw)-.5)
  freq=mid*pro%aw%d+pro%aw%o
  nz(iw)=calc_work(freq/6.285714,min_zsamp,max_zsamp,cycles,samples,dz(:,iw))
  mywork(iw)=io+nz(iw)
end do

lowest=10000.

#ifdef SEP_MPI
    if(pro%mpi_ithread==0) then
#endif


!Find the plan that equalizes the work as much as possible
do itry=1,10*nw
	cost_node=0;
	
  not_distrib=.true.

	do while(any(not_distrib))
    inode=minloc(cost_node)
    call random_number(myrand)
    igroup=nint(myrand*nw)
    igroup=min(nw,max(1,igroup))
    if(not_distrib(igroup)) then
      cost_node(inode(1))=cost_node(inode(1))+mywork(igroup)
      not_distrib(igroup)=.false.
      cur_plan(igroup)=inode(1)
    end if
  end do
 
  mycost=(maxval(cost_node)-minval(cost_node))/minval(cost_node)
  if(mycost<lowest) then
    lowest=mycost
    low_plan=cur_plan 
  end if
end do


if(pro%verb) write(pro%io_unit,*) "cost differential",lowest

#ifdef SEP_MPI
end if
call MPI_Bcast(low_plan,size(low_plan),MPI_INTEGER, 0, MPI_COMM_WORLD,ierr)
#endif
low_plan=low_plan-1
pro%iw_thread=low_plan

!Figure out how many groups local iw is holding
pro%nw_local=0
do i=1,size(low_plan)
  if(low_plan(i)==pro%mpi_ithread) pro%nw_local=pro%nw_local+1
end do

!Allocate the number of zs
allocate(pro%z_samp(pro%nw_local))

j=0
allocate(pro%iw_first(pro%nw_local))



  if(rite_sample) then
    call to_history("n1",maxval(nz),"sample_z")
    call to_history("n2",size(low_plan),"sample_z")
  end if

do i=1,size(low_plan)
  if(low_plan(i)==pro%mpi_ithread) then
    j=j+1
    pro%iw_first(j)=(do_freq_block(i)-1)*pro%nws+1
		allocate(pro%z_samp(j)%dz(nz(i)))
    pro%z_samp(j)%dz(1:nz(i))=dz(1:nz(i),i)
  end if
  if(rite_sample) then
    call srite("sample_z",dz(:,i),nz(i)*4)
    do k=nz(i)+1,maxval(nz)
      call srite("sample_z",0.,4)
    end do
  end if
end do

end subroutine












  !----------------------------------------------------------------
end module wei_work




