7
7
# nor does it submit to any jurisdiction.
8
8
9
9
from typing import Any
10
- from typing import Dict
11
- from typing import List
12
- from typing import Optional
13
10
from typing import Tuple
14
11
from typing import Union
15
12
16
13
import numpy as np
17
14
from numpy .typing import NDArray
18
15
19
16
from earthkit .meteo import constants
20
- from earthkit . meteo . thermo import specific_gas_consant
17
+
21
18
22
19
def pressure_at_model_levels (
23
20
A : NDArray [Any ], B : NDArray [Any ], sp : Union [float , NDArray [Any ]]
@@ -82,22 +79,29 @@ def pressure_at_model_levels(
82
79
# calculate alpha
83
80
alpha = np .zeros (new_shape_full )
84
81
85
- alpha [1 :, ...] = 1.0 - p_half_level [1 :- 1 , ...] / (p_half_level [2 :, ...] - p_half_level [1 :- 1 , ...]) * delta [1 :, ...]
82
+ alpha [1 :, ...] = (
83
+ 1.0 - p_half_level [1 :- 1 , ...] / (p_half_level [2 :, ...] - p_half_level [1 :- 1 , ...]) * delta [1 :, ...]
84
+ )
86
85
87
86
# pressure at highest half level <= 0.1
88
87
if np .any (p_half_level [0 , ...] <= PRESSURE_TOA ):
89
88
alpha [0 , ...] = 1.0 # ARPEGE choice, ECMWF IFS uses log(2)
90
89
# pressure at highest half level > 0.1
91
90
else :
92
- alpha [0 , ...] = 1.0 - p_half_level [0 , ...] / (p_half_level [1 , ...] - p_half_level [0 , ...]) * delta [0 , ...]
91
+ alpha [0 , ...] = (
92
+ 1.0 - p_half_level [0 , ...] / (p_half_level [1 , ...] - p_half_level [0 , ...]) * delta [0 , ...]
93
+ )
93
94
94
95
# calculate pressure on model full levels
95
96
# TODO: is there a faster way to calculate the averages?
96
97
# TODO: introduce option to calculate full levels in more complicated way
97
- p_full_level = np .apply_along_axis (lambda m : np .convolve (m , np .ones (2 ) / 2 , mode = "valid" ), axis = 0 , arr = p_half_level )
98
+ p_full_level = np .apply_along_axis (
99
+ lambda m : np .convolve (m , np .ones (2 ) / 2 , mode = "valid" ), axis = 0 , arr = p_half_level
100
+ )
98
101
99
102
return p_full_level , p_half_level , delta , alpha
100
103
104
+
101
105
def relative_geopotential_thickness (alpha : NDArray [Any ], q : NDArray [Any ], T : NDArray [Any ]) -> NDArray [Any ]:
102
106
"""Calculates the geopotential thickness w.r.t the surface on model full-levels.
103
107
@@ -115,13 +119,15 @@ def relative_geopotential_thickness(alpha: NDArray[Any], q: NDArray[Any], T: NDA
115
119
ndarray
116
120
geopotential thickness of model full-levels w.r.t. the surface
117
121
"""
122
+ from earthkit .meteo import specific_gas_constant
118
123
119
- R = calc_specific_gas_constant (q )
124
+ R = specific_gas_constant (q )
120
125
dphi = np .cumsum (np .flip (alpha * R * T , axis = 0 ), axis = 0 )
121
126
dphi = np .flip (dphi , axis = 0 )
122
127
123
128
return dphi
124
129
130
+
125
131
def pressure_at_height_level (
126
132
height : float , q : NDArray [Any ], T : NDArray [Any ], sp : NDArray [Any ], A : NDArray [Any ], B : NDArray [Any ]
127
133
) -> Union [float , NDArray [Any ]]:
@@ -154,7 +160,7 @@ def pressure_at_height_level(
154
160
tdphi = height * constants .g
155
161
156
162
# pressure(-related) variables
157
- p_full , p_half , _ , alpha = model_level_pressure (A , B , sp )
163
+ p_full , p_half , _ , alpha = pressure_at_model_levels (A , B , sp )
158
164
159
165
# relative geopot. thickness of full levels
160
166
dphi = relative_geopotential_thickness (alpha , q , T )
0 commit comments