program dgauss

use sep
use igrad1

implicit none
logical                       :: lobes
integer			:: i1,n1,n2,stat
real                          :: o1,o2,d1,d2
real, allocatable	          :: in(:,:),out(:)

call sep_init()
call sep_begin_prog()

call from_history("n1",n1); call from_history("d1",d1); call from_history("o1",o1)
call from_history("n3",n2); call from_history("d3",d2); call from_history("o3",o2)
call to_history("n1",n1); call to_history("d1",d1); call to_history("o1",o1)
call to_history("n2",1); call to_history("d2",1.); call to_history("o2",0.)
call to_history("n3",n2); call to_history("d3",d2); call to_history("o3",o2)
call to_history("eout",4);call to_history("label3","")
call from_param("lobes",lobes,.true.)
allocate(in(n1,n2),out(n1))

call sreed("in",in,4*n1*n2)
do i1=1,n2
   out=0
   stat = igrad1_lop ( .false., .false., in(:,i1), out)
   in(:,i1)=0
   stat = igrad1_lop (  .true., .false., in(:,i1), out)
end do
if (not(lobes)) then
   if ( abs(minval(in(1:int(n1/2.),int(n2/2.)))) > abs(maxval(in(1:int(n1/2.),int(n2/2.))))) then  
      do i1=1,n2
         where(in(1:int(n1/2.),i1) > 0) in(1:int(n1/2.),i1) =0.
      end do
   else
      do i1=1,n2
         where(in(1:int(n1/2.),i1) < 0) in(1:int(n1/2.),i1) =0.
      end do
   end if
   if ( abs(minval(in(int(n1/2.):n1,int(n2/2.)))) < abs(maxval(in(int(n1/2.):n1,int(n2/2.))))) then  
      do i1=1,n2
         where(in(int(n1/2.)+1:n1,i1) < 0) in(int(n1/2.):n1,i1) =0.
      end do
   else
      do i1=1,n2
         where(in(int(n1/2.)+1:n1,i1) > 0) in(int(n1/2.):n1,i1) =0.
      end do
   end if

end if

call srite("out",in,4*n1*n2)

deallocate(in,out)
call sep_end_prog()
end program
