module brutermo
use sep
use brute_general
implicit none
real, allocatable, private       :: scan(:,:,:,:), dso(:)
integer, private                 :: nodes
contains
!==================================================================
subroutine objfcn(x)
!
integer                :: i(1)
REAL                   :: x(ns1,ns3),rhocorr(ns1,ns3,nfunc),rho0(ns1,ns3),y1,y2,fac
REAL,allocatable       :: f(:,:),f1(:,:),f2(:,:)
integer                :: i1,i3,i4,inode,i5
integer                :: min1,min3,max1,max3,im1,im3
integer, external      :: omp_get_thread_num,omp_get_num_threads

 if (corridor) then
    call sreed("rho0",rho0,4*ns1*ns3)
    allocate(f(nfunc,nodes),f1(nfunc,nodes),f2(nfunc,nodes))
    do i4=1,nfunc
       do i3=1,ns3
          do i1=1,ns1
             y1=-perc1+2*perc1*(i4-1)/(nfunc-1)
             y2=-perc2+2*perc2*(i4-1)/(nfunc-1)
             fac=y1+(y2-y1)/(ns1-1)*(i1-1)
             rhocorr(i1,i3,i4)=rho0(i1,i3)*(1+fac/100)
          end do
       end do
    end do
    do i3=1,ns3
       do i1=1,ns1
          where(rhocorr(i1,i3,:)>(ns4-1)*ds4+os4)
             rhocorr(i1,i3,:)=rho0(i1,i3)  
          end where
          where(rhocorr(i1,i3,:)<os4)
             rhocorr(i1,i3,:)=rho0(i1,i3)  
          end where
          rhocorr(i1,i3,:)=int((rhocorr(i1,i3,:)-os4)/ds4)+1
          if(mod(i1,20).eq.1.and.mod(i3,5).eq.1) write(0,*) i1,i3,minval(rhocorr(i1,i3,:)),maxval(rhocorr(i1,i3,:))
       end do
    end do
    
 else
    allocate(f(ns4,nodes),f1(ns4,nodes),f2(ns4,nodes))
 end if
!$OMP parallel do private(inode,i3,min3,max3,im3,i1,min1,max1,im1,i4,i,i5) schedule(dynamic,1) num_threads(nodes)
do i3=1,ns3
   inode = omp_get_thread_num()+1
   min3=max(1,i3-jx)
   max3=min(ns3,i3+jx)
   f(:,inode)=0;f1(:,inode)=0;f2(:,inode)=0;
   do i1=1,ns1
      min1=max(1,i1-jz)
      max1=min(ns1,i1+jz)
      do im3=min3,max3
         do im1=min1,max1
            if (corridor) then
               do i4=1,nfunc
                  i5=int(rhocorr(i1,i3,i4))
!                  f1(i4,inode) = f1(i4,inode) + eps1*sum((scan(im1,:,im3,i5)*dso)**2)/ns2
!                  f2(i4,inode) = f2(i4,inode) - eps2*(scan(im1,int(ns2/2.)+1,im3,i5))**2
                  f1(i4,inode) = f1(i4,inode) + eps1*sum(abs(scan(im1,:,im3,i5)*dso))/(ns2)
                  f2(i4,inode) = f2(i4,inode) - eps2*abs(scan(im1,int(ns2/2.)+1,im3,i5))
                  f(i4,inode) = f(i4,inode) + f1(i4,inode) + f2(i4,inode)
               end do
            else
               do i4=1,ns4
                  f1(i4,inode) = f1(i4,inode) + eps1*sum((scan(im1,:,im3,i4)*dso)**2)/ns2
                  f2(i4,inode) = f2(i4,inode) - eps2*(scan(im1,int(ns2/2.)+1,im3,i4))**2
!                  f1(i4,inode) = f1(i4,inode) + eps1*sum(abs(scan(im1,:,im3,i4)*dso))/(ns2)
!                  f2(i4,inode) = f2(i4,inode) - eps2*abs(scan(im1,int(ns2/2.)+1,im3,i4))
                  f(i4,inode) = f(i4,inode) + f1(i4,inode) + f2(i4,inode)
               end do
            end if
         end do
      end do
      i=minloc(f(:,inode))
      if (corridor) then
         x(i1,i3) = (rhocorr(i1,i3,i(1))-1)*ds4+os4
      else
         x(i1,i3) = (i(1)-1)*ds4+os4
      end if
   end do
   if (mod(i3,50).eq.1) write(0,*) "Picked rmo for cmp ",i3,' from ', ns3
end do
!$OMP END PARALLEL DO
deallocate(f,f1,f2)
RETURN
end subroutine objfcn
!==================================================================
subroutine kernel(x)
real     :: x(:),xcan(ns1,ns2,ns3),rho
integer  :: i1,i2,i3,i4
do i1=1,ns1;do i3=1,ns3
  rho=x((i3-1)*ns1+i1)
  i4=(rho-os4)/ds4+1.5
  xcan(i1,:,i3)=scan(i1,:,i3,i4)!*dso
end do;end do
call to_history("n1",ns1,"xscan");call to_history("d1",ds1,"xscan");call to_history("o1",os1,"xscan");
call to_history("n2",ns2,"xscan");call to_history("d2",ds2,"xscan");call to_history("o2",os2,"xscan");
call to_history("n3",ns3,"xscan");call to_history("d3",ds3,"xscan");call to_history("o3",os3,"xscan");
call srite("xscan",xcan,4*ns1*ns2*ns3)
end subroutine 
!==================================================================
subroutine sa_init()
  integer           :: i
  integer, external :: omp_get_thread_num,omp_get_num_threads

if (node.eq.9999) then
  !$OMP PARALLEL
  nodes = omp_get_num_threads()
  !$OMP END PARALLEL
else
  nodes=node
end if
write(3,*) "****************"
write(3,*) "* NUM.THREADS = ",nodes
write(3,*) "****************"

  allocate(scan(ns1,ns2,ns3,ns4),dso(ns2))
  call sreed("scan",scan,4*ns1*ns2*ns3*ns4)
  do i=1,ns2
     dso(i)=(i-1)*ds2+os2
  end do
  dso=abs(dso)

end subroutine
!==================================================================
subroutine sa_close()
  deallocate(scan,dso)
end subroutine
!==================================================================
end module
