# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module tri_smooth_3d_mod
  use triangle
    implicit none
  interface triangle_smooth
    module procedure triangle_smooth_2d
    module procedure triangle_smooth_3d
  end interface
  contains
  logical function triangle_smooth_2d(field,nbox)
    real, pointer :: stor1(:)
    real :: field(:,:)
    integer :: nbox(2)
    integer :: i1,i2,i3
    triangle_smooth_2d=.false.
    allocate(stor1(size(field,1)))
    call triangle_init(nbox(1),size(field,1))
    do i2=1,size(field,2)
      stor1=field(:,i2)
      i1=triangle_lop(.false.,.false.,stor1,field(:,i2))
    end do
    deallocate(stor1)
    call triangle_close()
    allocate(stor1(size(field,2)))
    call triangle_init(nbox(2),size(field,2))
    do i1=1,size(field,1)
      stor1=field(i1,:)
      i2=triangle_lop(.false.,.false.,stor1,field(i1,:))
    end do
    deallocate(stor1)
    call triangle_close()
    triangle_smooth_2d=.true.
  end function
  logical function triangle_smooth_3d(field,nbox)
    real, pointer :: stor1(:)
    real :: field(:,:,:)
    integer :: nbox(3)
    integer :: i1,i2,i3
    triangle_smooth_3d=.false.
    do i3=1,size(field,3)
      if (.not. triangle_smooth_2d(field(:,:,i3),nbox(1:2))) then
        write(0,*) "        ERROR: Smooothing triangle smooth 2d"
        return
      end if
    end do
    if (size(field,3).ne.1) then
      allocate(stor1(size(field,3)))
      call triangle_init(nbox(3),size(field,3))
      do i2=1,size(field,2)
        do i1=1,size(field,1)
          stor1=field(i1,i2,:)
          i3=triangle_lop(.false.,.false.,stor1,field(i1,i2,:))
        end do
      end do
      deallocate(stor1)
      call triangle_close()
    end if
    triangle_smooth_3d=.true.
  end function
end module
