Числовая точность расчетов климатологии в Python

В настоящее время я пытаюсь воссоздать результаты статьи (https://www.researchgate.net/publication/309723672_Evidence_for_wave_resonance_as_a_key_mechanism_for_generating_high-amplitude_quasi-stationary_waves_in_boreal_summer) для моей магистерской диссертации. Я вычисляю меридиональное (градусы северной широты) распределение квадрата меридионального волнового числа волны Россби (определенный вид течения в атмосфере) для ряда дней. Это значение зависит только от средних зональных (градусы западного и восточного) ветра и их первой и второй меридиональной производной, а также от зонального волнового числа двумерной волны Россби. Волну такого типа можно представить, например, как волну на барабане, только в сферической среде, в атмосфере. Я использую python 3.6.5 и подозреваю, что проблема заключается в числовой точности, однако я не уверен.

Я читал другие потоки, касающиеся числовой точности, и наткнулся, например, на этот: точность . Однако я еще не пробовал этого, потому что стараюсь не писать свои собственные тригонометрические функции. Кроме того, поскольку мне приходится обрабатывать довольно большой объем данных, я стараюсь не замедлять свой код. Из экспериментов я обнаружил, что математическая библиотека не более точна, чем библиотека numpy, в отношении тригонометрических функций.

Вот фрагмент кода, который меня беспокоит:

Lat = np.linspace(0,90,37)
MeridWN = np.zeros((29,36), dtype='float64')        
        ######################################################
        #define Meridional wavenumber, l^2
        for i in range(5,28,10):
            for j in range(36):
                MeridWN[i,j] = (((2*EarthRot*np.cos(Lat[j]*np.pi/180.0)**3.0)/(EarthRad*ZonMeanZonWiNH[j]))-
                   ((np.cos(Lat[j]*np.pi/180.)**2.)/(EarthRad**2.0*ZonMeanZonWiNH[j]))*
                   ZonMeanZonWiMeridGradGrad[j]+
                   ((np.sin(Lat[j]*np.pi/180.)*np.cos(Lat[j]*np.pi/180.))/(EarthRad**2.0*ZonMeanZonWiNH[j]))*
                   ZonMeanZonWiMeridGrad[j]+(1./(EarthRad**2.0))-(ZonWN[i]/EarthRad)**2.0)
                MeridWNMerge[i,x,j] = MeridWN[i,j]

Индекс i относится к диапазону зональных волновых чисел, x — день (этот фрагмент взят из более крупного цикла, работающего в течение нескольких дней), а j — положение широты. Для вычисления производных я использую функции градиента numpy следующим образом:

ZonMeanZonWiMeridGrad = np.gradient(ZonMeanZonWiNH,np.linspace(0,90,37))
ZonMeanZonWiMeridGradGrad = np.gradient(ZonMeanZonWiMeridGrad,np.linspace(0,90,37))

Это формула для расчета квадрата меридионального волнового числа (l), где Omega - вращение Земли, Phi - широтное положение, a - радиус Земли, U - среднее зональное, усредненное по зонам, а k - зональное волновое число, в моем случае массив в диапазоне от 5,5 до 8,5.

Это сравнение моего среднего зонального поля ветра (внизу) с июня по август и один на бумаге (вверху), что указывает на то, что у нас одинаковые данные, и это не проблема. Цветовые шкалы немного различаются, однако наиболее характерные особенности профиля ветра очень похожи, и небольшие различия не должны давать такой разницы профиль меридионального волнового числа (для k = 7), где рисунок из бумаги снова сверху, а мой снизу. Здесь цветовые шкалы снова разные, но, тем не менее, следует уловить большое структурное сходство. Как видите, я имею дело с очень маленькими числами, что привело к подозрению в том, что у меня неточный в числовом отношении код.
Если хотите, я могу загрузить весь свой код, однако я думаю, что для обсуждения точности этого достаточно. Я пытаюсь решить эту проблему около 2 недель, пробуя все разные изменения в моем коде, некоторые из них были хорошими изменениями, но ни одно из них не дало желаемого результата.

Заранее спасибо,
Томас