#!/usr/bin/env python3
# scan_bndy_cnst.py — varre BNDY.file.* e CNST.file de uma rodada do Eta
# procurando valores nao-fisicos (NaN/Inf) e checando a estrutura vertical.
#
# Criado no diagnostico do crash do topo 1 hPa (2026-08-01): TABLE.f com
# pl=pt gravava Inf/NaN nas tabelas BM (PTBL/TTBL/qs0/sqs/sthe/the0) do
# CNST — regs ~122 (2332 B) e ~132 (82136 B) p/ itb=76/jtb=134.
# Apos regerar o CNST com o TABLE.f corrigido, esses registros devem sair
# com NaN/Inf=0.
#
# Uso:
#   python3 scan_bndy_cnst.py <dir_rodada> [IM JM LM [NSOIL]]
#   (defaults: IM=477 JM=719 LM=50 NSOIL=8; espera <dir>/CNST.file e
#    um <dir>/BNDY.file.* — o primeiro que achar)
import glob
import os
import struct
import sys

import numpy as np

run = sys.argv[1] if len(sys.argv) > 1 else "."
IM = int(sys.argv[2]) if len(sys.argv) > 2 else 477
JM = int(sys.argv[3]) if len(sys.argv) > 3 else 719
LM = int(sys.argv[4]) if len(sys.argv) > 4 else 50
NSOIL = int(sys.argv[5]) if len(sys.argv) > 5 else 8
LB = 2 * IM + JM - 3

def recs(path):
    f = open(path, "rb")
    while True:
        h = f.read(4)
        if len(h) < 4:
            break
        n = struct.unpack("<i", h)[0]
        d = f.read(n)
        f.read(4)
        yield d

def stats(name, a):
    nan = np.isnan(a).sum()
    inf = np.isinf(a).sum()
    print(f"{name:14s} min={np.nanmin(a):14.6e} max={np.nanmax(a):14.6e} NaN={nan} Inf={inf}")
    return nan + inf

bad_total = 0

bndys = sorted(glob.glob(os.path.join(run, "BNDY.file.*")))
if bndys:
    print(f"========= {bndys[0]} =========")
    g = recs(bndys[0])
    hdr = next(g)
    print("header run,idat,ihrst,tboco:", struct.unpack("<i3iif", hdr[:24]))
    blk = 0
    while True:
        try:
            bchr = np.frombuffer(next(g), dtype="<f4")
        except StopIteration:
            break
        print(f"--- bloco {blk}: BCHR={bchr}")
        pdb = np.frombuffer(next(g), dtype="<f4").reshape(2, LB)
        bad_total += stats("PDB valor", pdb[0])
        bad_total += stats("PDB tend", pdb[1])
        for nm in ("TB", "QB", "UB", "VB", "Q2B", "CWMB"):
            d = np.frombuffer(next(g), dtype="<f4").reshape(2, LM, LB)
            bad_total += stats(nm + " valor", d[0])
            bad_total += stats(nm + " tend", d[1])
        blk += 1
else:
    print("(sem BNDY.file.* em", run, ")")

cnst = os.path.join(run, "CNST.file")
print(f"\n========= {cnst} =========")
sizes, data = [], []
for d in recs(cnst):
    sizes.append(len(d))
    data.append(d)
print("total regs:", len(sizes))

r1 = data[0]
print("reg1 bytes:", len(r1), " JM deduzido:", (len(r1) + 200) / 48)
lmh = np.frombuffer(data[1], dtype="<i4").reshape(JM, IM)
lmv = np.frombuffer(data[2], dtype="<i4").reshape(JM, IM)
print(f"LMH min/max: {lmh.min()}/{lmh.max()}  LMV min/max: {lmv.min()}/{lmv.max()}")
for k, nm in enumerate(("HBM2", "VBM2", "VBM3", "SM", "SICE")):
    bad_total += stats(nm, np.frombuffer(data[3 + k], dtype="<f4").reshape(JM, IM))
htm_bad = sum(((np.frombuffer(data[8 + l], dtype="<f4") != 0)
               & (np.frombuffer(data[8 + l], dtype="<f4") != 1)).sum() for l in range(LM))
vtm_bad = sum(((np.frombuffer(data[8 + LM + l], dtype="<f4") != 0)
               & (np.frombuffer(data[8 + LM + l], dtype="<f4") != 1)).sum() for l in range(LM))
print("HTM valores != 0/1:", htm_bad, "  VTM != 0/1:", vtm_bad)

rc = data[8 + 2 * LM]
scal = struct.unpack("<10f", rc[:40])
print("dy,cpgfv,en,ent,r,pt,tddamp,f4d,f4q,ef4t =", ["%g" % s for s in scal])
off = 40
detac = np.frombuffer(rc, dtype="<f4", count=LM, offset=off); off += 4 * LM
rdeta = np.frombuffer(rc, dtype="<f4", count=LM, offset=off); off += 4 * LM
aeta  = np.frombuffer(rc, dtype="<f4", count=LM, offset=off); off += 4 * LM
off += 4 * LM  # f4q2
eta   = np.frombuffer(rc, dtype="<f4", count=LM + 1, offset=off); off += 4 * (LM + 1)
dfl   = np.frombuffer(rc, dtype="<f4", count=LM + 1, offset=off)
print("PT no arquivo:", scal[5])
print("DETAC: min=%g max=%g zeros=%d" % (detac.min(), detac.max(), (detac == 0).sum()))
print("RDETA: min=%g max=%g NaN/Inf=%d" % (rdeta.min(), rdeta.max(), (~np.isfinite(rdeta)).sum()))
print("AETA mono?", bool(np.all(np.diff(aeta) > 0)), " ETA[0],ETA[-1]=", eta[0], eta[-1])
print("DFL topo..sfc:", dfl[0], "...", dfl[-1], " mono decresc?", bool(np.all(np.diff(dfl) < 0)))

print("\n--- varredura NaN/Inf de TODOS os registros do CNST ---")
achou = False
for k, n in enumerate(sizes):
    a = np.frombuffer(data[k], dtype="<f4")
    bad = int((~np.isfinite(a)).sum())
    if bad:
        achou = True
        bad_total += bad
        print(f"reg {k:3d} bytes={n:9d}  NaN/Inf={bad}  min={np.nanmin(a):.4g} max={np.nanmax(a):.4g}")
if not achou:
    print("nenhum registro com NaN/Inf  ✓  (tabelas BM sadias)")

print("\ntotal geral de valores nao-finitos:", bad_total)
sys.exit(1 if bad_total else 0)