#!/usr/bin/env python3
"""Gera deta_70_1mb — LM=70, topo PT=1 hPa, BASE = deta_50_1mb (PBL identica).

Desenho (2026-08-07): compromisso 50L->70L que gasta os 20 niveis extras SO nos
pontos sensiveis de altos niveis, mantendo a PBL do 50L (13 camadas do
deta_50_25mb VERBATIM em dp) => acoplamento de superficie identico ao 50L
(comparacoes 50x70 isolam a vertical superior).
  - Zona S (L1-11, 1->44 hPa): dlnp 0.60->0.21; 5 camadas inteiras acima de
    10 hPa (50L: 4); topo-buffer.
  - Zona TTL (L12-22, 44->144 hPa): bolsao da inversao, dlnp min 0.083
    centrado no cold point (81-96 hPa; dz 0.48-0.51 km; 50L: ~0.85 km).
  - Zona M (L23-57, 144->888 hPa): jato 145-310 com dp 17-22 (50L: 26-36);
    pico 26.3 hPa em ~487 (50L: 42) — qualidade proxima do deta_60.
  - Zona B (L58-70): 13 camadas de PBL do 50L VERBATIM.
Juntas 0.98/1.04/0.94. Requer PARMETA LM=70 + co2.70_1mb (gen_co2) + CNST.
Custo ~ +40%% vs LM=50.
"""
import numpy as np, struct, os

SUP = '/home/jorge/Eta_support_data/static/fix/deta_files'
raw = open(os.path.join(SUP, 'deta_50_25mb'), 'rb').read()
nrec = struct.unpack('<i', raw[:4])[0]
d50  = np.frombuffer(raw[4:204], dtype='<f4').astype(float)
t50  = np.frombuffer(raw[204:4+nrec], dtype='<i4')          # LDUM(51)
dp50 = d50*(1000.-25.)

PT = 1.; PS = 1000.; COL = PS-PT
dlnpS = np.array([0.60,0.51,0.44,0.38,0.34,0.31,0.28,0.26,0.24,0.22,0.21])
dlnpT = np.array([0.17,0.14,0.115,0.095,0.086,0.083,0.083,0.086,0.095,0.107,0.118])
dlnp  = np.concatenate([dlnpS, dlnpT])
pS  = PT*np.exp(np.concatenate([[0], np.cumsum(dlnp)]))
dpS = np.diff(pS)
nbot = 13
dpB  = dp50[-nbot:]
nmid = 70-len(dpS)-nbot
span = (PS-dpB.sum())-pS[-1]
base = np.linspace(dpS[-1]*1.04, dpB[0]/0.94, nmid)
x = np.linspace(0, 1, nmid)
shape = np.sin(np.pi*x)**1.5
dpM = base + (span-base.sum())/shape.sum()*shape

dp   = np.concatenate([dpS, dpM, dpB])
deta = (dp/COL).astype(np.float32)          # residuo ~1e-7: interp.f forca eta(LM+1)=1.
ldum = np.round(np.interp(np.linspace(0, 1, 71),
                          np.linspace(0, 1, len(t50)), t50)).astype('<i4')

rec = deta.tobytes() + ldum.tobytes()
out = os.path.join(SUP, 'deta_70_1mb')
open(out, 'wb').write(struct.pack('<i', len(rec)) + rec + struct.pack('<i', len(rec)))

pint = np.concatenate([[PT], PT+np.cumsum(dp)])
print(f"gravado {out}: LM=70 PT=1hPa  sum={deta.sum(dtype=np.float32):.7f}")
print(f"{'L':>2} {'p_top':>8} {'p_bot':>8} {'dp_hPa':>7} {'DETA':>10}")
for L in range(70):
    print(f"{L+1:>2} {pint[L]:8.2f} {pint[L+1]:8.2f} {dp[L]:7.2f} {deta[L]:10.6f}")