Issue with CAPE Calculation

[Stone-Shaw]

I am encountering a problem where the calculated CAPE value is very large, but the TlnP plot shows it is actually quite small.
Version: 1.6.3
What I have tried:

1. Checked the MetPy source code and found different results when running step-by-step versus calling the function directly.
2. Removed the data above 100 hPa found different results again, but it is still different from step-by-step result.

I am unsure what the issue might be.

from metpy.units import units
from metpy.calc import parcel_profile 
from metpy.calc import lfc,el
from metpy.calc.thermo import _find_append_zero_crossings
from metpy.calc import dewpoint_from_relative_humidity
from metpy.calc import cape_cin
from metpy.calc.tools import _less_or_close, _greater_or_close
from metpy import constants as mpconsts
import numpy as np

P = [1000.,  975.,  950.,  925.,  900.,  875.,  850.,  825.,  800.,
        775.,  750.,  700.,  650.,  600.,  550.,  500.,  450.,  400.,
        350.,  300.,  250.,  225.,  200.,  175.,  150.,  125.,  100.,
         70.,   50.,   30.,   20.,   10.,    7.,    5.,    3.,    2.,
          1.] 
T = [ 26.601715 ,  25.321442 ,  23.64029  ,  22.3544   ,  21.61026  ,
        20.53189  ,  19.3378   ,  18.081207 ,  16.891754 ,  15.679352 ,
        14.605133 ,  11.797516 ,   8.922485 ,   4.584381 ,   0.8770447,
        -2.8042297,  -7.650421 , -13.071838 , -19.521637 , -27.249115 ,
       -36.776184 , -41.891373 , -47.43611  , -55.09149  , -63.31491  ,
       -70.46591  , -74.436005 , -73.86389  , -62.693848 , -52.720825 ,
       -47.223602 , -41.97264  , -34.355484 , -29.62732  , -15.841644 ,
        -4.863159 , -11.00827  ]

R = [8.7531734e-01, 8.7178373e-01, 9.2426902e-01, 9.7321486e-01,
       9.7158509e-01, 9.6182829e-01, 9.7018987e-01, 9.7185606e-01,
       9.5916349e-01, 9.4788438e-01, 9.2271966e-01, 9.3404025e-01,
       8.3645624e-01, 8.2847106e-01, 8.3991808e-01, 8.0008334e-01,
       7.7837539e-01, 7.8936952e-01, 8.8870806e-01, 8.9821351e-01,
       9.3429244e-01, 1.0480621e+00, 9.7989351e-01, 9.9125290e-01,
       1.0505068e+00, 9.9932271e-01, 6.5298730e-01, 2.0224983e-01,
       2.8747728e-02, 4.7798837e-03, 1.6737868e-03, 4.7243937e-04,
       1.4527883e-04, 6.6689492e-05, 1.0917375e-05, 2.7924066e-06,
       2.4962396e-06]

P_value_unit = P * units.hPa
T_value_unit = T * units.degC
R_value_unit = R * units.dimensionless
Td_value_unit = dewpoint_from_relative_humidity(T_value_unit, R_value_unit)

def cape_cin_steps(pressure,temperature,dewpoint, parcel_profile):
	lfc_pressure, _ = lfc(P_value_unit, T_value_unit, Td_value_unit,
						parcel_temperature_profile=parcel_profile, which='bottom')

	# If there is no LFC, no need to proceed.
	if np.isnan(lfc_pressure):
		print("LFC_p is None")
		cape = units.Quantity(0, 'J/kg')
		cin  = units.Quantity(0, 'J/kg')
	else:
		el_pressure, _ = el(pressure, temperature, dewpoint,
						parcel_temperature_profile=parcel_profile, which='top')

		if np.isnan(el_pressure):
				el_pressure = pressure[-1]
		else:
				el_pressure = el_pressure

		y = (parcel_profile - temperature).to(units.degK)
		x, y = _find_append_zero_crossings(np.copy(pressure), y)

		# CAPE
		# Only use data between the LFC and EL for calculation
		p_mask = _less_or_close(x, lfc_pressure) & _greater_or_close(x, el_pressure)
		x_clipped = x[p_mask].magnitude
		y_clipped = y[p_mask].magnitude
		cape = (mpconsts.Rd
				* units.Quantity(np.trapz(y_clipped, np.log(x_clipped)), 'K')).to(units('J/kg'))
		# CIN
		# Only use data between the surface and LFC for calculation
		p_mask = _greater_or_close(x, lfc_pressure)
		x_clipped = x[p_mask].magnitude
		y_clipped = y[p_mask].magnitude
		cin = (mpconsts.Rd
			* units.Quantity(np.trapz(y_clipped, np.log(x_clipped)), 'K')).to(units('J/kg'))
		if cin > units.Quantity(0, 'J/kg'):
			cin = units.Quantity(0, 'J/kg')
	return(cape,cin)

parcel_profile_value_unit = parcel_profile(P_value_unit, T_value_unit[0], Td_value_unit[0]).to('degC')
cape,cin = cape_cin_steps(P_value_unit, T_value_unit, Td_value_unit, parcel_profile_value_unit)
print(f"My CAPE:{cape:.2f}	MY CIN:{cin:.2f}")

#* metpy CAPE
CAPE,CIN = cape_cin(P_value_unit, T_value_unit, Td_value_unit, parcel_profile_value_unit)
print(f"Metpy CAPE:{CAPE:.2f}	 CIN:{CIN:.2f}")

output:

My CAPE:864.72 joule / kilogram	MY CIN:-12.96 joule / kilogram
Metpy CAPE:-114565.61 joule / kilogram	 CIN:-8.88 joule / kilogram

when I Removed the data above 100 hPa:

index = -10
P = P[:index]
T = T[:index]
R = R[:index]

output:

My CAPE:864.72 joule / kilogram	MY CIN:-12.96 joule / kilogram
Metpy CAPE:913.12 joule / kilogram	 CIN:-8.88 joule / kilogram