!===============================================================================
!  shapiro_gef_v3 : filtro de Shapiro 2D (ordem 1..4) para campos 2D do pós Eta
!
!  Uso:  ./shapiro_v3 [ordem] [npass] [arq_in] [irec] [arq_out]
!        defaults: 2 500 fort.10 1 mslp_filt_v3.bin
!        irec = 0  -> filtra TODOS os registros do arquivo (saída com o mesmo
!                     layout; basta trocar o DSET no .ctl original)
!
!  Diferenças em relação à v2:
!   - undef é preenchido por dilatação (média dos vizinhos válidos, iterada)
!     antes da filtragem e restaurado na saída -> a ordem cheia vale em todo o
!     domínio, sem degradação de resposta junto ao quadro de undef
!   - borda do arranjo por reflexão (halo de 4 pontos)
!   - irec=0 processa todos os registros
!  Mantido da v2: sinal correto para toda ordem, sem periodicidade, recl
!  portável, dupla precisão na acumulação, tabela de resposta.
!
!  Dimensionamento para as ondas de ~7 dx da PSLM 1 km (vetor de onda zonal):
!     ordem 2 x 500 passadas: R(7dx)~0 R(10dx)=0.01 R(15dx)=0.39 R(20dx)=0.74 R(30dx)=0.94
!     (2 x 150 zera 7dx mas deixa ondulação residual de ~10 dx)
!===============================================================================
program shapiro_v3
  implicit none
  integer, parameter :: nx = 906, ny = 441, nh = 4
  real(4),  parameter :: undef_in  = 1.e20, undef_out = -9999.
  integer  :: order = 2, npass = 500, irec = 1
  character(len=256) :: fin = 'fort.10', fout = 'mslp_filt_v3.bin', arg
  real(4), allocatable :: fio(:,:)
  real(8), allocatable :: f(:,:), g(:,:)
  logical, allocatable :: ok(:,:)
  real(8) :: coef(0:8, 0:4), pi, resp
  integer :: k, n, p, reclen, nrec, r, r1, r2, nfill
  integer(8) :: fsize

  if (command_argument_count() >= 1) then; call get_command_argument(1, arg); read(arg,*) order; end if
  if (command_argument_count() >= 2) then; call get_command_argument(2, arg); read(arg,*) npass; end if
  if (command_argument_count() >= 3) call get_command_argument(3, fin)
  if (command_argument_count() >= 4) then; call get_command_argument(4, arg); read(arg,*) irec;  end if
  if (command_argument_count() >= 5) call get_command_argument(5, fout)
  if (order < 1 .or. order > nh) stop 'ordem suportada: 1..4'

  coef = 0.d0
  do n = 1, 4
    do k = 0, 2*n
      coef(k, n) = (-1.d0)**k * binom(2*n, k) * (-1.d0)**n / 4.d0**n
    end do
  end do

  pi = 4.d0*atan(1.d0)
  print '(a,i2,a,i5,a)', 'Shapiro ordem ', order, ', ', npass, ' passadas. Resposta por comprimento de onda:'
  do k = 4, 30, 2
    resp = (1.d0 - sin(pi/real(k,8))**(2*order))**npass
    print '(a,i3,a,f8.4)', '   L = ', k, ' dx  ->  R = ', resp
  end do

  allocate(fio(nx,ny), f(1-nh:nx+nh, 1-nh:ny+nh), g(1-nh:nx+nh, 1-nh:ny+nh), ok(nx,ny))
  inquire(iolength=reclen) fio
  inquire(file=trim(fin), size=fsize)
  nrec = int(fsize / (int(nx,8)*int(ny,8)*4_8))
  if (irec == 0) then; r1 = 1; r2 = nrec; else; r1 = irec; r2 = irec; end if
  print '(a,i4,a,i4,a,i4)', 'registros no arquivo: ', nrec, '   processando ', r1, ' a ', r2

  open(10, file=trim(fin),  access='direct', form='unformatted', recl=reclen, status='old',     action='read')
  open(11, file=trim(fout), access='direct', form='unformatted', recl=reclen, status='replace', action='write')

  do r = r1, r2
    read(10, rec=r) fio
    ok = abs(fio) < 0.5*undef_in .and. fio /= undef_out
    f = 0.d0
    f(1:nx,1:ny) = fio
    call fill_undef(nfill)
    do p = 1, npass
      call reflect(f);  call pass_x(f, g)
      call reflect(g);  call pass_y(g, f)
    end do
    fio = real(f(1:nx,1:ny), 4)
    where (.not. ok) fio = undef_out
    write(11, rec=r-r1+1) fio
    print '(a,i3,a,i7,a,i6,a,2f12.4)', 'rec ', r, ': validos ', count(ok), '  iter.fill ', nfill, &
          '  min/max saida ', minval(fio, ok), maxval(fio, ok)
  end do
  close(10); close(11)
  print *, 'escrito ', trim(fout)

contains

  real(8) function binom(nn, kk)
    integer, intent(in) :: nn, kk
    integer :: m
    binom = 1.d0
    do m = 1, kk
      binom = binom * real(nn - kk + m, 8) / real(m, 8)
    end do
  end function binom

  ! preenche undef por dilatação iterativa (média dos vizinhos válidos, 8 vizinhos)
  subroutine fill_undef(niter)
    integer, intent(out) :: niter
    logical :: done(nx,ny), newp(nx,ny)
    real(8) :: s
    integer :: i, j, di, dj, c
    done = ok
    niter = 0
    do
      if (all(done)) exit
      niter = niter + 1
      if (niter > nx+ny) stop 'fill_undef: nao convergiu'
      newp = .false.
      do j = 1, ny
        do i = 1, nx
          if (done(i,j)) cycle
          s = 0.d0; c = 0
          do dj = -1, 1
            do di = -1, 1
              if (i+di < 1 .or. i+di > nx .or. j+dj < 1 .or. j+dj > ny) cycle
              if (done(i+di,j+dj)) then; s = s + f(i+di,j+dj); c = c + 1; end if
            end do
          end do
          if (c > 0) then; g(i,j) = s / c; newp(i,j) = .true.; end if
        end do
      end do
      if (.not. any(newp)) stop 'fill_undef: campo sem pontos validos'
      where (newp) f(1:nx,1:ny) = g(1:nx,1:ny)
      done = done .or. newp
    end do
  end subroutine fill_undef

  ! halo por reflexão: a(0)=a(2), a(-1)=a(3), ...; a(nx+1)=a(nx-1), ...
  subroutine reflect(a)
    real(8), intent(inout) :: a(1-nh:nx+nh, 1-nh:ny+nh)
    integer :: h
    do h = 1, nh
      a(1-h, 1:ny)  = a(1+h, 1:ny)
      a(nx+h, 1:ny) = a(nx-h, 1:ny)
    end do
    do h = 1, nh
      a(:, 1-h)  = a(:, 1+h)
      a(:, ny+h) = a(:, ny-h)
    end do
  end subroutine reflect

  subroutine pass_x(a, b)
    real(8), intent(in)  :: a(1-nh:nx+nh, 1-nh:ny+nh)
    real(8), intent(out) :: b(1-nh:nx+nh, 1-nh:ny+nh)
    integer :: i, j, k
    real(8) :: s
    do j = 1, ny
      do i = 1, nx
        s = 0.d0
        do k = 0, 2*order
          s = s + coef(k,order) * a(i-order+k, j)
        end do
        b(i,j) = a(i,j) - s
      end do
    end do
  end subroutine pass_x

  subroutine pass_y(a, b)
    real(8), intent(in)  :: a(1-nh:nx+nh, 1-nh:ny+nh)
    real(8), intent(out) :: b(1-nh:nx+nh, 1-nh:ny+nh)
    integer :: i, j, k
    real(8) :: s
    do j = 1, ny
      do i = 1, nx
        s = 0.d0
        do k = 0, 2*order
          s = s + coef(k,order) * a(i, j-order+k)
        end do
        b(i,j) = a(i,j) - s
      end do
    end do
  end subroutine pass_y

end program shapiro_v3