# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module semb_mod
  use sep
  use sep3d_struct_mod
  implicit none
  real, pointer,private,save :: numer(:),denom(:),count(:)
  integer,private,save :: nwind1
  integer, private, external :: sep_thread_num
  contains
  logical function init_semb(input)
    type(sep3d) :: input
    call from_param("nsmooth",nwind1,min(5,input%n(1)))
    write(0,*) "INI"
    allocate(numer(input%n(1)),denom(input%n(1)),count(input%n(1)))
    init_semb=.true.
  end function
  logical function semb_clean()
    deallocate(numer,denom,count)
    semb_clean=.true.
  end function
  logical function semb_it(buffer,semb)
    complex :: buffer(:,:,:,:,:)
    real :: semb(:,:,:)
    real :: numsum,densum,amp
    integer :: i4,i5,i1,i2,i3
    semb_it=.false.
    semb=0.
    do i5=1,size(buffer,5)
      do i4=1,size(buffer,4)
        numer=0
        denom=0
        count=0.
        do i3=1,size(buffer,3)
          do i2=1,size(buffer,2)
            do i1=1,size(buffer,1)
              amp=real(buffer(i1,i2,i3,i4,i5))
              if (amp .ne.0.) then
                numer(i1)= numer(i1)+amp
                denom(i1)= denom(i1)+amp**2
                count(i1)= count(i1)+1
              end if
            end do
          end do
        end do
        do i1=1,size(numer)
          numsum=0
          densum=0
          do i2=max(1,i1-nwind1),min(i1+nwind1,size(numer))
            if (count(i2) > 0.) then
              numsum= numsum+numer(i2)**2
              densum= densum+denom(i2)*count(i2)
            end if
          end do
          if (densum>0.) then
            semb(i1,i4,i5)=numsum/densum
          end if
        end do
      end do
    end do
    semb_it=.true.
  end function
end module
