Numerical instability: pi() returns NaN for physically reasonable atmospheric profile

Summary

The pi() function returns NaN values for a physically reasonable atmospheric profile extracted from ERA5 reanalysis data. The same algorithm succeeds with very similar profiles from neighboring dates/conditions.

Environment

tcpyPI version: 1.4
Python version: 3.12.11
OS: macOS (Darwin)

Reproduction

The following minimal script reproduces the issue:

"""Minimal reproduction of tcpyPI NaNs issue.

We construct two atmospheric profiles - one that causes tcpyPI to return NaNs, and
a control profile that works successfully. They're associated with subsequent days
in the same location (so surface conditions and pressure levels are nearly identical).
We then call tcpyPI.pi directly with these profiles.
"""

import numpy as np
from tcpyPI.pi import pi

# Extracted the exact problematic atmospheric profile from our data
# This is from dayofyear=287, year=1995 at lat=34.75°N, lon=311.5°E

# Surface conditions (converted to tcpyPI expected units)
SSTC = 25.267  # Sea surface temp in Celsius
MSL = 1012.9  # Mean sea level pressure in hPa


# (pressure in hPa (decreasing), temperature in Celsius, specific humidity in g/kg)
# tcpyPI expects decreasing order for pressure levels
profile_problematic = np.array(
    [
        [1000.0, 20.266, 9.65],
        [975.0, 18.11, 9.517],
        [950.0, 15.978, 9.459],
        [925.0, 13.94, 9.205],
        [900.0, 12.576, 7.784],
        [875.0, 11.213, 6.561],
        [850.0, 9.669, 5.952],
        [825.0, 8.084, 5.517],
        [800.0, 6.516, 4.976],
        [775.0, 4.993, 4.332],
        [750.0, 4.374, 3.102],
        [700.0, 4.93, 1.207],
        [650.0, 1.852, 1.521],
        [600.0, -2.047, 1.263],
        [550.0, -6.917, 0.857],
        [500.0, -13.015, 0.556],
        [450.0, -19.698, 0.33],
        [400.0, -26.881, 0.211],
        [350.0, -34.178, 0.104],
        [300.0, -42.822, 0.073],
        [250.0, -51.573, 0.037],
        [225.0, -54.68, 0.027],
        [200.0, -57.207, 0.018],
        [175.0, -58.943, 0.01],
        [150.0, -60.751, 0.005],
        [125.0, -63.255, 0.004],
        [100.0, -63.277, 0.003],
        [70.0, -61.5, 0.003],
        [50.0, -58.453, 0.003],
        [30.0, -53.327, 0.003],
        [20.0, -50.156, 0.003],
        [10.0, -46.602, 0.003],
        [7.0, -42.283, 0.003],
        [5.0, -35.389, 0.003],
        [3.0, -22.349, 0.003],
        [2.0, -14.944, 0.004],
        [1.0, -11.886, 0.004],
    ]
)

P = profile_problematic[:, 0]
TC = profile_problematic[:, 1]
R = profile_problematic[:, 2]

# tcpyPI options
pi_options = dict(CKCD=0.9, ascent_flag=0, diss_flag=1, ptop=50, miss_handle=1)

print()
print("=" * 50)
print("TEST: Testing a profile that causes numerical instability...")
print("\t(dayofyear=287, year=1995)")
print(
    f"\tSurface conditions: SST: {SSTC:.3f} °C, MSL: {MSL:.1f} hPa, "
    f"pressure levels from {P[0]} to {P[-1]} hPa"
)

# Call tcpyPI potential intensity function directly
VMAX, PMIN, IFL, TO, OTL = pi(
    SSTC,  # Sea surface temperature (°C)
    MSL,  # Mean sea level pressure (hPa)
    P,  # Pressure levels (hPa)
    TC,  # Temperature profile (°C)
    R,  # Specific humidity profile (g/kg)
    **pi_options,
)

print("\t...")
print("\tResults:")
print(f"\tVMAX (max wind speed): {VMAX:.3f} m/s")
print(f"\tPMIN (min pressure):   {PMIN:.3f} hPa")
print(f"\tIFL (flag):            {IFL}")
print(f"\tTO (outflow temp):     {TO:.3f} K")
print(f"\tOTL (outflow level):   {OTL:.3f} hPa")

# Check for NaNs
if np.isnan(VMAX):
    print("\tREPRODUCTION SUCCESSFUL!")
    print("\ttcpyPI returned NaN for VMAX with seemingly valid input data")
else:
    print("\tUNEXPECTED: tcpyPI succeeded with this profile")

# Working profile data (P is the same as above, temperature and humidity differ)
profile_working = np.array(
    [
        [1000.0, 21.326, 11.215],
        [975.0, 19.192, 11.011],
        [950.0, 17.315, 10.556],
        [925.0, 15.676, 9.931],
        [900.0, 14.286, 8.639],
        [875.0, 12.957, 7.58],
        [850.0, 11.495, 6.97],
        [825.0, 10.019, 6.439],
        [800.0, 8.575, 5.734],
        [775.0, 6.993, 5.198],
        [750.0, 5.373, 4.752],
        [700.0, 3.13, 3.122],
        [650.0, 0.534, 2.562],
        [600.0, -2.573, 1.442],
        [550.0, -6.616, 0.865],
        [500.0, -12.233, 0.61],
        [450.0, -18.954, 0.42],
        [400.0, -26.098, 0.292],
        [350.0, -33.337, 0.184],
        [300.0, -40.639, 0.112],
        [250.0, -49.373, 0.057],
        [225.0, -52.236, 0.03],
        [200.0, -55.706, 0.017],
        [175.0, -57.323, 0.01],
        [150.0, -60.084, 0.005],
        [125.0, -65.224, 0.004],
        [100.0, -65.651, 0.003],
        [70.0, -63.457, 0.003],
        [50.0, -59.522, 0.003],
        [30.0, -53.504, 0.003],
        [20.0, -49.926, 0.003],
        [10.0, -45.925, 0.003],
        [7.0, -41.476, 0.003],
        [5.0, -34.68, 0.003],
        [3.0, -23.495, 0.003],
        [2.0, -16.575, 0.004],
        [1.0, -12.898, 0.004],
    ]
)

TC_working = profile_working[:, 1]
R_working = profile_working[:, 2]
SSTC_working = 25.292
MSL_working = 1011.4

print()
print("CONTROL TEST: Testing a known working profile...")
print("\t(dayofyear=286, year=1995)")
print(
    f"\tSurface conditions: SST: {SSTC_working:.3f} °C, MSL: {MSL_working:.1f} hPa, "
    f"pressure levels from {P[0]} to {P[-1]} hPa"
)


VMAX_work, PMIN_work, IFL_work, TO_work, OTL_work = pi(
    SSTC_working, MSL_working, P, TC_working, R_working, **pi_options
)

print("\t...")
print("\tResults:")
print(f"\tVMAX (max wind speed): {VMAX_work:.3f} m/s")
print(f"\tPMIN (min pressure):   {PMIN_work:.3f} hPa")
print(f"\tIFL (flag):            {IFL_work}")
print(f"\tTO (outflow temp):     {TO_work:.3f} K")
print(f"\tOTL (outflow level):   {OTL_work:.3f} hPa")

if np.isnan(VMAX_work):
    print("\tControl test also failed - broader issue")
else:
    print(f"\tControl test succeeded: VMAX={VMAX_work:.1f} m/s")

Expected vs Actual Behavior

Expected: The algorithm should return valid potential intensity values or a clear error message.

Actual: Returns NaN for all output variables with IFL=0, indicating internal algorithm failure.

Additional Context

The atmospheric profile comes from ERA5 reanalysis data (October 14, 1995)
All input values are within physically reasonable ranges
Very similar profiles from neighboring dates work successfully
This appears to be a numerical instability rather than invalid input data
The profile has slightly cooler temperatures and lower humidity compared to working cases, but still within normal meteorological bounds

This issue was discovered during climate attribution analysis where 1 out of 930 historical atmospheric profiles caused regression failures due to NaN PI values.

Might be the same as #77, where I isolated something very similar. I suspect it's a problem with the fixed point iteration having to do with the discretization of atmospheric levels resulting in 2-cycles (flip-flopping between two values) rather than 1-cycles (i.e. convergence). I'm tracking this internally as AT-990.

@maresb

I'be begun working on a branch to trace this problem and believe I have tracked it down to this line:

tcpyPI/src/tcpyPI/pi.py

Line 601 in 3c1d4ef

while (np.abs(PNEW-PMOLD) > 0.5):

With the block in pi.py:

# loop until convergence or bail out
    while (np.abs(PNEW-PMOLD) > 0.5):
        
        #
        #   ***  Find CAPE at radius of maximum winds   ***
        #
        TP=T[NK]
        PP=min([PM,1000.0])
        # find the mixing ratio with the average of the lowest level pressure and MSL
        RP=constants.EPS*R[NK]*MSL/(PP*(constants.EPS+R[NK])-R[NK]*MSL)
        result = cape(TP,RP,PP,T,R,P,ascent_flag,ptop,miss_handle)
        CAPEM = result[0]
        IFLAG = result[3]
        # if the CAPE function tripped a different flag, set the output IFL to it
        if (IFLAG != 1):
            IFL=int(IFLAG)
        
        #
        #  ***  Find saturation CAPE at radius of maximum winds    ***
        #  *** Note that TO and OTL are found with this assumption ***
        #
        TP=SSTK
        PP=min([PM,1000.0])
        RP=utilities.rv(ES0,PP)
        result = cape(TP,RP,PP,T,R,P,ascent_flag,ptop,miss_handle)
        CAPEMS, TOMS, LNBS, IFLAG = result
        # if the CAPE function tripped a flag, set the output IFL to it
        if (IFLAG != 1):
            IFL=int(IFLAG)
        # Store the outflow temperature and level of neutral bouyancy at the outflow level (OTL)
        TO=TOMS   
        OTL=LNBS
        # Calculate the proxy for TC efficiency (BE02, EQN. 1-3)
        RAT=SSTK/TO
        # If dissipative heating is "off", TC efficiency proxy is set to 1.0 (BE02, pg. 3)
        if (diss_flag == 0):
            RAT=1.0
        
        #
        #  ***  Initial estimate of pressure at the radius of maximum winds  ***
        #
        RS0=RP
        # Lowest level and Sea-surface Density Temperature (E94, EQN. 4.3.1 and 6.3.7)
        TV0=utilities.Trho(T[NK],R[NK],R[NK])
        TVSST=utilities.Trho(SSTK,RS0,RS0)
        # Average Surface Density Temperature, e.g. 1/2*[Tv(Tsfc)+Tv(sst)]
        TVAV=0.5*(TV0+TVSST)
        # Converge toward CAPE*-CAPEM (BE02, EQN 3-4)
        CAT=(CAPEM-CAPEA)+0.5*CKCD*RAT*(CAPEMS-CAPEM)
        CAT=max([CAT,0.0])
        # Iterate on pressure
        PNEW=MSL*np.exp(-CAT/(constants.RD*TVAV))
        
        #
        #   ***  Test for convergence (setup for possible next while iteration)  ***
        #
        # store the previous step's pressure       
        PMOLD=PM
        # store the current step's pressure
        PM=PNEW
        # increase iteration count in the loop
        NP += 1
        
        #
        #   ***   If the routine does not converge, set IFL=0 and return missing PI   ***
        #
        if (NP > 200)  or (PM < 400):
            VMAX=np.nan
            PMIN=np.nan
            IFL=0
            TO=np.nan
            OTL=np.nan
            return(VMAX,PMIN,IFL,TO,OTL)
    
    # Once converged, set potential intensity at the radius of maximum winds
    CATFAC=0.5*(1.+1/constants.b)
    CAT=(CAPEM-CAPEA)+CKCD*RAT*CATFAC*(CAPEMS-CAPEM)
    CAT=max([CAT,0.0])

This doesn't converge within 200 iterations of NP (and note that from the flag there's just a generic failure flag here, i.e IFL=0). The PM values are oscillating with PNEW-PMOLD just greater than 0.5 hPa.

What does an error of 0.5hPa mean? I've probed this a little bit and it's about the same as a 0.5 m/s error---so it's a 1-to-1 scaling between hPa and m/s errors over range of reasonable inputs. In other words, our VMAX precision is only as good as the precision on this convergence algorithm.

The first change I propose to tcpyPI is simply to update the flag such that IFL=2 when it doesn't converge, to match the other failure to converge flag. I think this change should be made regardless of what else we do.

To address the convergence issue, we could either:

A) keep outputting as NaNs, recognizing this lack of convergence happens sometimes and we aren't willing to accept an error > 0.5 m/s. This is a reasonable scientific approach. In this case, we would need to deal with the NaNs in our attribution pipeline somehow.

B) update the algorithm so that if the PNEW-PMOLD difference is > 0.5 hPa but < 1.0 hPa, we take a final PM of (PNEW-PMOLD)/2, and then move the code along. This would minimize the errors between them and allow a valid VMAX output, and would only affect the profiles where this convergence is failing (and no where else).

C) change the convergence threshold to a wider range (0.5 --> 1.0 hPa)... I worry about this option because I don't know how it will affect ALL the PI values we regularly calculate with this algorithm, if we uniformly weaken their precision.

@maresb @alansammarone please let me know what you think.

@dgilford

That's great - thanks for looking into this @dgilford !

I don't have a good sense of the tradeoffs here, specially between B and C. From the POV of the pipeline, we can probably make A work as well, at least up to some threshold like we talked about. I'll leave it to you and Ben since I don't have anything useful to say heheh but I'll keep you posted about my own investigation

@dgilford

Hi @dgilford, I believe the problem, at least for #77, is here:

tcpyPI/src/tcpyPI/pi.py

Lines 243 to 256 in 3c1d4ef

    
               INB=0 
        
               for j in range(nlvl-1, jmin, -1): 
        
                   if (TVRDIF[j] > 0): 
        
                       INB=max([INB,j]) 
        
               # CHECK: Is the LNB higher than the surface? If not, return zero CAPE   
        
               if (INB==0): 
        
                   CAPED=0 
        
                   TOB=T[0] 
        
                   LNB=P[INB] 
        
           #         TOB=np.nan 
        
                   LNB=0 
        
                   # Return the unconverged values 
        
                   return(CAPED,TOB,LNB,IFLAG)

This causes a discontinuity in CAPEM vs PP around where the parcel first becomes buoyant, so there is no root to be found.

I don't know enough about the physics to suggest a good fix.

The data in the following plot is from #77.

To bridge the gap on this for now, we have decided to move forward with option A, that is:

A) keep outputting as NaNs, recognizing this lack of convergence happens sometimes and we aren't willing to accept an error > 0.5 m/s.

As noted above,

The first change I propose to tcpyPI is simply to update the flag such that IFL=2 when it doesn't converge, to match the other failure to converge flag. I think this change should be made regardless of what else we do.

So I'm proposing that I will go ahead and make the above change, i.e. set IFL=2, and merge into main.

For now the convergence limit in L601 will remain hard-coded at +/- 0.5 hPa.

#81 has merged a fix that makes verbose IFL=2 when the pressure code fails to converge, addressing this particular issue for the above example profile. Otherwise, the code is now (still) working as intended, returning NaN when the profile conditions are such that the pressure algorithm fails to converge.

I'm closing this for now, with the convergence limit remaining hard-coded at +/- 0.5 hPa, which still leads to NaN for sometimes physically "reasonable" (but non-convergent) profiles.

Navigation Menu

Repository navigation

Numerical instability: pi() returns NaN for physically reasonable atmospheric profile #80

Summary

Environment

Reproduction

Expected vs Actual Behavior

Additional Context

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

	INB=0
	for j in range(nlvl-1, jmin, -1):
	if (TVRDIF[j] > 0):
	INB=max([INB,j])

	# CHECK: Is the LNB higher than the surface? If not, return zero CAPE
	if (INB==0):
	CAPED=0
	TOB=T[0]
	LNB=P[INB]
	# TOB=np.nan
	LNB=0
	# Return the unconverged values
	return(CAPED,TOB,LNB,IFLAG)

Numerical instability: pi() returns NaN for physically reasonable atmospheric profile #80

Description

Summary

Environment

Reproduction

Expected vs Actual Behavior

Additional Context

Activity

maresb commented on Sep 12, 2025

dgilford commented on Sep 17, 2025

alansammarone commented on Sep 17, 2025

maresb commented on Sep 18, 2025

dgilford commented on Sep 24, 2025

dgilford commented on Oct 7, 2025

Metadata

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Issue actions