# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module gsmooth_mod
  implicit none
  contains
  subroutine gauss_smooth(model,nwind1,nwind2,nwind3,sd)
    integer n1,n2,i1,i2,n3,i3
    real, allocatable :: input(:,:,:),output(:,:,:),gauss(:,:,:)
    integer :: nwind1,nwind2,nwind3,c1,c2,c3
    integer :: bd3,ed3,bd2,ed2,bd1,ed1
    integer :: bg3,eg3,bg2,eg2,bg1,eg1
    integer :: l1,l2,l3
    real :: sm,d1,d3,sd,d2,model(:,:,:)
    n1=size(model,1)
    n2=size(model,2)
    n3=size(model,3)
    allocate(input(n1,n2,n3),output(n1,n2,n3))
    allocate(gauss(nwind1,nwind2,nwind3))
    input=model
    if (mod(nwind1,2).ne.1 .or. mod(nwind2,2).ne.1 .or. mod(nwind3,2&
      &).ne.1) then
      call erexit("nwind must be odd")
    end if
    c1=ceiling(nwind1/2.)
    c2=ceiling(nwind2/2.)
    c3=ceiling(nwind3/2.)
    d1=-log(.0000001)
    d2=((c1-1)**2+(c2-1)**2+(c3-1)**2)**.5
    sd= (d2/d1**.5)
    do i3=1,nwind3
      d3=(c3-i3)**2
      do i2=1,nwind2
        d2=(c2-i2)**2
        do i1=1,nwind1
          d1=((c1-i1)**2+d2+d3)**.5
          gauss(i1,i2,i3)=exp(-((d1)/sd)**2)
        end do
      end do
    end do
    l1=nwind1-c1
    l2=nwind2-c2
    l3=nwind3-c3
    output=0.
!gauss=1
    do i3=1,n3
      bd3=max(1,i3-l3)
      ed3=min(n3,i3+l3)
      bg3=c3-(i3-bd3)
      eg3=c3+(ed3-i3)
      do i2=1,n2
        bd2=max(1,i2-l2)
        ed2=min(n2,i2+l2)
        bg2=c2-(i2-bd2)
        eg2=c2+(ed2-i2)
        do i1=1,n1
          bd1=max(1,i1-l1)
          ed1=min(n1,i1+l1)
          bg1=c1-(i1-bd1)
          eg1=c1+(ed1-i1)
          sm=sum(gauss(bg1:eg1,bg2:eg2,bg3:eg3))
          output(i1,i2,i3)=sum(input(bd1:ed1,bd2:ed2,bd3:ed3)*gauss&
            &(bg1:eg1,bg2:eg2,bg3:eg3))/sm
        end do
      end do
    end do
    model=output
    deallocate(input,output,gauss)
  end subroutine
end module
