program ab
 use quant_mod
 use sep
 implicit none
 integer :: n1,n2,i,k,j,ii,n3
 real ,allocatable :: buf1(:),buf2(:),temp(:)
 real :: pct1(100),pct2(100),md,f
 integer :: ic
 call sep_init()
 call from_history("n1",n1)
 call from_history("n2",n2)
 call from_history("n3",n3)
 allocate(buf1(n1*n2*n3),buf2(n1*n2*n3),temp(n1*n2*n3))
write(0,*) " aa",n1,n2,n3,size(buf1),size(buf2)

 call sreed("in",buf1,size(buf1)*4)
write(0,*) "DY1 "
 call sreed("other",buf2,size(buf2)*4)

write(0,*) " aa",n1,n2,n3
 temp=buf1
 do i=1,99
write(0,*) "PERCENTILE 1",i
   k=real(i)/100.*n1*n2*n3
   call quantinternal(k,n1*n2*n3,temp)
   pct1(i)=temp(k)
 end do

 
 temp=buf2
 do i=1,100
write(0,*) "PERCENTILE 2",i
   k=real(i)/100.*n1*n2*n3
   call quantinternal(k,n1*n2*n3,temp)
   pct2(i)=temp(k)
 end do

write(0,*) "WHAT THE",size(pct1),size(pct2)
call to_history("n1",99,"map")
call to_history("n2",2,"map")
call srite("map",pct1,396)
call srite("map",pct2,396)
write(0,*) "WHA2 THE",size(pct1),size(pct2)

 ii=n1*n2*n3/100
 do i=1,n1*n2*n3
if(mod(i,ii)==0) write(0,*) "finished",nint(i*100./(n1*n2*n3*1.))
   if(pct1(1)>buf1(i)) then
     buf1(i)=0.
   else if(buf1(i) > pct1(99))  then
     buf1(i)=100.
   else
     j=2
     do while(pct1(j)  < buf1(i)) 
       j=j+1
     end do
     f=(buf1(i)-pct1(j-1))/1.
     buf1(i)=f+j-1
   end if
 end do

 
 do i=1,n1*n2*n3
   if(pct2(1)>buf2(i)) then
     buf2(i)=0.
   else if(buf2(i) > pct2(99))  then
     buf2(i)=100.
   else
     j=2
     do while(pct2(j)  < buf2(i)) 
       j=j+1
     end do
     f=(buf2(i)-pct2(j-1))/1.
     buf2(i)=f+j-1
   end if
 end do
 buf1=buf1-buf2
 call srite("out", buf1,size(buf1)*4)

 
 




end program
