남극 빙하 빙붕 저면 융해 열염순환(Thermohaline) 3방정식 & 1,000Hz 0-Alloc C# 해수면 상승 질량 결손 역해석 엔진

using System;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Stem.EarthScience.Glaciology
{
    /// 

    /// 남극 빙붕 저면 융해 열염순환 3방정식 및 질량 결손 역해석 구조체 (0-Alloc C#)
    /// 
    public ref struct IceShelfBasalMeltEngine
    {
        private const double WaterDensity = 1028.0;      // kg/m^3 (해수 밀도)
        private const double IceDensity = 917.0;         // kg/m^3 (빙하 밀도)
        private const double WaterHeatCapacity = 3974.0; // J/(kg*K)
        private const double LatentHeatMelt = 3.34e5;    // J/kg
        private const double DragCoeffCd = 0.0025;       // 항력 계수
        private const double GammaT_factor = 0.011;      // 열 교환 감마 상수
        private const double GammaS_factor = 0.00031;    // 염분 교환 감마 상수
        
        // 빙점 강하 상수 (UNESCO 공식 단순화)
        private const double A = -0.0573; // °C / PSU
        private const double B = 0.0832;  // °C
        private const double C = -7.53e-8; // °C / Pa

        private readonly double _iceConductivityTerm; // W/m^2 (내부 빙하 전도 열류)
        private double _currentMeltRateMPerYear;
        private double _cumulativeIceLossTons;

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public IceShelfBasalMeltEngine(double internalIceTemperature, double iceThickness)
        {
            // 내부 얼음 전도 열류 근사: rho_i * c_i * kappa_i * (dT/dz)
            double tempGrad = (0.0 - internalIceTemperature) / Math.Max(iceThickness, 10.0);
            _iceConductivityTerm = 917.0 * 2009.0 * 1.54e-6 * tempGrad;
            _currentMeltRateMPerYear = 0.0;
            _cumulativeIceLossTons = 0.0;
        }

        /// 

        /// 원격 심해 관측 수온, 염분, 해류 유속, 수심 압력을 기반으로 3방정식 뉴턴-랩슨 반복해를 수렴 (0-Alloc)
        /// 
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public void SolveBoundaryStep(
            double farFieldTemp, 
            double farFieldSalinity, 
            double currentVelocity, 
            double depthMeters, 
            double areaSquareMeters,
            double dtSeconds,
            out double meltRateMPerSec, 
            out double boundaryTemp, 
            out double boundarySalinity)
        {
            double pressurePa = WaterDensity * 9.80665 * depthMeters;
            double uStar = Math.Sqrt(DragCoeffCd) * Math.Max(currentVelocity, 0.01);
            double gammaT = uStar * GammaT_factor;
            double gammaS = uStar * GammaS_factor;

            // 3방정식 경계 염분 Sb 초기 추정
            double sb = farFieldSalinity;
            double tb = (A * sb) + B + (C * pressurePa);
            double meltRate = 0.0;

            // 빠른 수렴 뉴턴-랩슨 4회 루프 (0-Alloc 고정)
            for (int iter = 0; iter < 4; iter++)
            {
                tb = (A * sb) + B + (C * pressurePa);
                double qThermal = WaterDensity * WaterHeatCapacity * gammaT * (farFieldTemp - tb);
                meltRate = Math.Max(0.0, (qThermal - _iceConductivityTerm) / (IceDensity * LatentHeatMelt));
                
                // 염분 수지식 Sb 갱신: gammaS * (Sw - Sb) = (rho_i / rho_w) * meltRate * Sb
                double denom = gammaS + (IceDensity / WaterDensity) * meltRate;
                if (denom > 1e-9)
                {
                    sb = (gammaS * farFieldSalinity) / denom;
                }
            }

            meltRateMPerSec = meltRate;
            boundaryTemp = tb;
            boundarySalinity = sb;
            
            // 연간 환산 융해율 (m/yr) 및 누적 손실량 (톤)
            _currentMeltRateMPerYear = meltRate * 365.25 * 86400.0;
            double stepMassLossKg = IceDensity * meltRate * areaSquareMeters * dtSeconds;
            _cumulativeIceLossTons += (stepMassLossKg / 1000.0);
        }

        public double CurrentMeltRateMPerYear => _currentMeltRateMPerYear;
        public double CumulativeIceLossTons => _cumulativeIceLossTons;
    }
}