#!/usr/bin/env python3
"""Gera deta_80_1mb — distribuicao vertical LM=80 para topo PT=1 hPa.

Desenho v2 (2026-08-06, mesma metodologia do deta_50_1mb + bolsao da inversao):
  - Zona S (L1-12, 1->47 hPa): Delta(ln p) 0.55->0.21, dz 4.2->1.3 km monotonico
    (sem degrau refletor p/ ondas); topo espesso = buffer da tampa (+CODAMP/DDAMP).
    5 camadas inteiras acima de 10 hPa.
  - Zona TTL (L13-24, 47->146 hPa): BOLSAO de refinamento na inversao de T em
    altos niveis (cold point tropical ~90-100 hPa): Delta(ln p) desce a 0.08 no
    centro -> dz 0.47-0.68 km em 55-135 hPa (9 camadas), entrada/saida graduais.
    Resolve o salto de N^2 da tropopausa (reflexao/duto de ondas, cirrus TTL).
  - Zona M (L25-56, 146->830 hPa): parametrica, juntas 1.04/0.95; jato 150-300
    com dp 18-21 hPa (dz 0.5-0.74, cobre a tropopausa extratropical do dominio);
    pico 24.6 hPa em ~470 hPa (= deta_60 operacional).
  - Zona B (L57-80, 830->1000 hPa): as 24 ultimas camadas do deta_60 VERBATIM
    em dp (PBL/superficie identica a operacional de 60 niveis; 1a camada ~19 m).
Requisitos de uso: PARMETA LM=80 (recompilar etafcst), tabelas CO2 co2.80_1mb,
CNST regerado; custo dinamica ~ +60%% vs LM=50.
Formato: 1 registro unformatted LE = DETA(80) R4 + LDUM(81) I4 (dummy;
interpolado do tail do deta_60 so por plausibilidade — interp.f le e ignora).
"""
import numpy as np, struct, os

SUP = '/home/jorge/Eta_support_data/static/fix/deta_files'
raw = open(os.path.join(SUP, 'deta_60'), 'rb').read()
nrec = struct.unpack('<i', raw[:4])[0]
d60  = np.frombuffer(raw[4:244], dtype='<f4').astype(float)
t60  = np.frombuffer(raw[244:4+nrec], dtype='<i4')          # LDUM(61)
PT60, PS = 25., 1000.
dp60 = d60*(PS-PT60)

PT = 1.; COL = PS-PT
dlnpS = np.array([0.55,0.46,0.40,0.35,0.31,0.28,0.26,0.24,0.23,0.22,0.21,0.21])
dlnpT = np.array([0.17,0.14,0.115,0.10,0.09,0.082,0.078,0.078,0.082,0.09,0.102,0.118])
dlnp  = np.concatenate([dlnpS, dlnpT])
pS  = PT*np.exp(np.concatenate([[0], np.cumsum(dlnp)]))
dpS = np.diff(pS)
nbot = 24
dpB  = dp60[-nbot:]
nmid = 80-len(dpS)-nbot
span = (PS-dpB.sum())-pS[-1]
base = np.linspace(dpS[-1]*1.04, dpB[0]/0.95, nmid)   # 18.3 -> 17.2: quase plana; o bump da o pico

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, 81),
                          np.linspace(0, 1, len(t60)), t60)).astype('<i4')

rec = deta.tobytes() + ldum.tobytes()
out = os.path.join(SUP, 'deta_80_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=80 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(80):
    print(f"{L+1:>2} {pint[L]:8.2f} {pint[L+1]:8.2f} {dp[L]:7.2f} {deta[L]:10.6f}")