! scan_init.f90 — varre o preproc.init do initbc procurando campos
! nao-fisicos (NaN, PD<=0/absurdo, T fora de faixa, Q<0/absurdo, ref<=0)
! que expliquem SIGSEGV na 1a chamada da radiacao (E290/FST88).
!
! Compilar (mesmo ambiente Intel do initbc):
!   ifort -O0 -traceback scan_init.f90 -o scan_init
! Uso:
!   ./scan_init <caminho>/preproc.init IM JM LM
program scan_init
  implicit none
  integer :: im, jm, lm, ios, n, i, j, l, nrep
  character(len=512) :: fname, a1, a2, a3
  logical(4) :: run
  integer(4) :: idat(3), ihrst, ntsd
  real(4), allocatable :: ueta(:,:,:), veta(:,:,:), teta(:,:,:), qeta(:,:,:)
  real(4), allocatable :: pd(:,:), phis(:,:), sm(:,:), ref(:,:)
  real(4), allocatable :: eta(:), deta(:), etal(:), zeta(:)
  real(4) :: pt

  call get_command_argument(1, fname)
  call get_command_argument(2, a1); read(a1,*) im
  call get_command_argument(3, a2); read(a2,*) jm
  call get_command_argument(4, a3); read(a3,*) lm

  allocate(ueta(im,jm,lm), veta(im,jm,lm), teta(im,jm,lm), qeta(im,jm,lm))
  allocate(pd(im,jm), phis(im,jm), sm(im,jm), ref(im,jm))
  allocate(eta(lm+1), deta(lm), etal(lm), zeta(lm+1))

  open(10, file=trim(fname), status='old', form='unformatted', iostat=ios)
  if (ios /= 0) stop 'erro abrindo arquivo'
  read(10, iostat=ios) run, idat, ihrst, ntsd, ueta, veta
  if (ios /= 0) stop 'erro lendo registro 1 (conferir IM/JM/LM)'
  read(10, iostat=ios) teta, qeta, pd, phis, sm, ref, eta, pt, deta, etal, zeta
  if (ios /= 0) stop 'erro lendo registro 2 (conferir IM/JM/LM)'
  close(10)

  write(*,'(a,3i5,a,i3,a,i6)') 'idat=', idat, '  ihrst=', ihrst, '  ntsd=', ntsd
  write(*,'(a,f10.2)') 'pt=', pt

  ! ---- PD: o unico campo sem salvaguarda no caminho lHydroInterp ----
  write(*,'(/a)') '=== PD (deve ser ~ps-pt: 20000..102000 Pa) ==='
  write(*,'(a,2es14.6)') 'min/max: ', minval(pd), maxval(pd)
  n = count(pd /= pd)
  write(*,'(a,i10)') 'NaN: ', n
  nrep = 0
  do j = 1, jm
    do i = 1, im
      if (pd(i,j) /= pd(i,j) .or. pd(i,j) <= 0. .or. pd(i,j) > 110000.) then
        nrep = nrep + 1
        if (nrep <= 20) write(*,'(a,2i6,es14.6,a,f10.1,a,f8.4)') &
          '  PD ruim em (i,j)=', i, j, pd(i,j), '  phis/g=', phis(i,j)/9.8, &
          '  ref=', ref(i,j)
      end if
    end do
  end do
  write(*,'(a,i10)') 'total PD fora de [0,110000]: ', nrep

  ! ---- TETA ----
  write(*,'(/a)') '=== TETA (radiacao exige [100,380) K) ==='
  write(*,'(a,2es14.6)') 'min/max: ', minval(teta), maxval(teta)
  write(*,'(a,i10)') 'NaN: ', count(teta /= teta)
  write(*,'(a,i10)') 'fora [100,380): ', count(teta < 100. .or. teta >= 380.)
  write(*,'(a,i10,a,i10)') 'no clamp 150: ', count(teta == 150.), &
    '   no clamp 325: ', count(teta == 325.)
  nrep = 0
  do l = 1, lm
    do j = 1, jm
      do i = 1, im
        if (teta(i,j,l) /= teta(i,j,l)) then
          nrep = nrep + 1
          if (nrep <= 10) write(*,'(a,3i6)') '  TETA NaN em (i,j,l)=', i, j, l
        end if
      end do
    end do
  end do

  ! ---- QETA ----
  write(*,'(/a)') '=== QETA (esperado [0, ~0.03] kg/kg) ==='
  write(*,'(a,2es14.6)') 'min/max: ', minval(qeta), maxval(qeta)
  write(*,'(a,i10)') 'NaN: ', count(qeta /= qeta)
  write(*,'(a,i10)') 'Q<0: ', count(qeta < 0.)
  write(*,'(a,i10)') 'Q>0.05: ', count(qeta > 0.05)

  ! ---- REF (o const.f faz res=1/res) ----
  write(*,'(/a)') '=== REF (0<ref<=1; const.f inverte) ==='
  write(*,'(a,2es14.6)') 'min/max: ', minval(ref), maxval(ref)
  write(*,'(a,i10)') 'ref<=0: ', count(ref <= 0.)

  ! ---- ventos (sanidade) ----
  write(*,'(/a)') '=== U/V (sanidade) ==='
  write(*,'(a,2es14.6)') 'U min/max: ', minval(ueta), maxval(ueta)
  write(*,'(a,2es14.6)') 'V min/max: ', minval(veta), maxval(veta)
  write(*,'(a,2i10)') 'U/V NaN: ', count(ueta /= ueta), count(veta /= veta)
end program scan_init