-
Notifications
You must be signed in to change notification settings - Fork 1
/
sth_model.py
66 lines (59 loc) · 2.42 KB
/
sth_model.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jun 8 14:27:58 2022
@author: [email protected]
"""
def sth_model(lat, doy, pvu=3.5):
'''
Sigmoid-model for tropopause height determination
INPUT
lat : latitude, in degrees [-90..90] (can be an array)
doy : day-of-year [1..366] (can be an array, same size as lat)
pvu : PVU value [1.5..3.5]
OUTPUT
z : tropopause height, in km
For details see:
Global empirical models for tropopause height determination
Pedro Mateus, Virgílio B. Mendes, and Carlos A.L. Pires
Instituto Dom Luiz, Faculdade de Ciências, Universidade de Lisboa, Lisbon, Portugal
Submitted to Atmospheric Research Journal
'''
from numpy import shape, array, polyfit, polyval, exp, pi, cos
if len( shape(lat) ) == 0:
lat = [lat]
doy = [doy]
if shape(lat) != shape(doy):
print('The lat and doy vector must be the same size!!!')
return
z = []
for i in range(0, len(lat)):
if lat[i] <= 0:
# SH
#a_0 a_1 a_2 a_3 a_4 a_5 pvu
a = array([\
[7.2086, 9.1185, -18.9363, 1.8065, 0.1067, -0.1859], # 1.5
[7.7749, 8.5516, -21.3633, 1.3565, 0.0879, -0.1640], # 2.0
[8.1797, 8.1455, -23.4839, 1.1464, 0.0798, -0.1491], # 2.5
[8.5055, 7.8175, -25.3146, 1.1133, 0.0819, -0.1370], # 3.0
[8.7796, 7.5417, -26.7300, 1.1474, 0.0870, -0.1255]]) # 3.5
else:
# NH
#a_0 a_1 a_2 a_3 a_4 a_5 pvu
a = array([\
[7.1454, 9.1847, 20.1173, -2.8115, 0.1509, -0.1863], # 1.5
[7.6501, 8.6773, 22.4300, -2.3154, 0.1334, -0.1642], # 2.0
[7.9925, 8.3329, 24.1731, -1.8069, 0.1082, -0.1493], # 2.5
[8.2649, 8.0584, 25.5622, -1.5508, 0.0950, -0.1372], # 3.0
[8.4990, 7.8229, 26.7284, -1.5804, 0.0980, -0.1257]]) # 3.5
a0, a1, a2, a3, a4, a5 = a[:,0], a[:,1], a[:,2], a[:,3], a[:,4], a[:,5]
d = (1+exp(-(lat[i]-a2)/a3))**a4
arg = 2*pi/365.25 * (doy[i]-28)
z_ = a0 + a1/d + a5 * cos(arg)
if pvu == 3.5:
z.append(z_[-1])
else:
p = polyfit([1.5, 2.0, 2.5, 3.0, 3.5], z_, 1)
z.append(polyval(p, pvu))
return array(z)