amescap.FV3_utils

FV3_utils contains internal Functions for processing data in MGCM output files such as vertical interpolation.

These functions can be used on their own outside of CAP if they are imported as a module:

from /u/path/FV3_utils import fms_press_calc

Third-party Requirements:

• numpy

• warnings

• scipy

Module Contents

Functions

MGStau_ls_lat(ls, lat)	Return the max altitude for the dust from "MGS scenario" from
MGSzmax_ls_lat(ls, lat)	Return the max altitude for the dust from "MGS scenario" from
UT_LTtxt(UT_sol[, lon_180, roundmin])	Returns the time in HH:MM:SS at a certain longitude.
add_cyclic(data, lon)	Add a cyclic (overlapping) point to a 2D array. Useful for azimuth
alt_KM(press[, scale_height_KM, reference_press])	Gives the approximate altitude [km] for a given pressure
area_meridional_cells_deg(lat_c, dlon, dlat[, ...])	Return area of invidual cells for a meridional band of thickness
area_weights_deg(var_shape, lat_c[, axis])	Return weights for averaging the variable.
areo_avg(VAR, areo, Ls_target, Ls_angle[, symmetric])	Return a value average over a central solar longitude
axis_interp(var_IN, x, xi, axis[, reverse_input, ...])	One dimensional linear/logarithmic interpolation along one axis.
azimuth2cart(LAT, LON, lat0[, lon0])	Azimuthal equidistant projection. Converts from latitude-longitude
broadcast(var_1D, shape_out, axis)	Broadcast a 1D array based on a variable's dimensions
cart_to_azimut_TR(u, v[, mode])	Convert cartesian coordinates or wind vectors to radians using azimuth angle.
compute_uneven_sigma(num_levels, N_scale_heights, ...)	Construct an initial array of sigma based on the number of levels
daily_to_average(varIN, dt_in[, nday, trim])	Bin a variable from an atmos_daily file format to the
daily_to_diurn(varIN, time_in)	Bin a variable from an atmos_daily file into the
dvar_dh(arr[, h])	Differentiate an array A[dim1, dim2, dim3...] w.r.t h. The
expand_index(Nindex, VAR_shape_axis_FIRST, axis_list)	Repeat interpolation indices along an axis.
find_n(X_IN, X_OUT[, reverse_input, modulo])	Maps the closest index from a 1D input array to a ND output array
find_n0(Lfull_IN, Llev_OUT[, reverse_input])	Return the index for the level(s) just below Llev_OUT.
fms_Z_calc(psfc, ak, bk, T[, topo, lev_type])	Returns the 3D altitude field [m] AGL (or above aeroid).
fms_press_calc(psfc, ak, bk[, lev_type])	Returns the 3D pressure field from the surface pressure and the
frontogenesis(U, V, theta, lon_deg, lat_deg[, R, spacing])	Compute the frontogenesis (local change in potential temperature
gauss_profile(x, alpha[, x0])	Return Gaussian line shape at x. This can be used to generate a
get_trend_2D(VAR, LON, LAT[, type_trend])	Extract spatial trends from the data. The output can be directly
interp_KDTree(var_IN, lat_IN, lon_IN, lat_OUT, lon_OUT)	Inverse distance-weighted interpolation using nearest neighboor for
layers_mid_point_to_boundary(pfull, sfc_val)	A general description for the layer boundaries is:
lin_interp(X_in, X_ref, Y_ref)	Simple linear interpolation with no dependance on scipy
lon180_to_360(lon)	Transform a float or an array from the -180/180 coordinate system
lon360_to_180(lon)	Transform a float or an array from the 0-360 coordinate system to
ls2sol(Ls_in)	Ls to sol converter.
mass_stream(v_avg, lat, level[, type, psfc, H, factor])	Compute the mass stream function:
mollweide2cart(LAT, LON)	Mollweide projection. Converts from latitude-longitude to
ortho2cart(LAT, LON, lat0[, lon0])	Orthographic projection. Converts from latitude-longitude to
polar2XYZ(lon, lat, alt[, Re])	Spherical to cartesian coordinate transformation.
polar_warming(T, lat[, outside_range])	Return the polar warming, following McDunn et al. 2013:
press_pa(alt_KM[, scale_height_KM, reference_press])	Gives the approximate altitude [km] for a given pressure
press_to_alt_atmosphere_Mars(Pi)	Return the altitude [m] as a function of pressure from the
ref_atmosphere_Mars_PTD(Zi)	Analytical atmospheric model for Martian pressure, temperature, and
regression_2D(X, Y, VAR[, order])	Linear and quadratic regression on the plane.
robin2cart(LAT, LON)	Robinson projection. Converts from latitude-longitude to cartesian
second_hhmmss(seconds[, lon_180])	Given the time [sec], return local true solar time at a
sfc_area_deg(lon1, lon2, lat1, lat2[, R])	Return the surface between two sets of latitudes/longitudes:
shiftgrid_180_to_360(lon, data)	This function shifts ND data from a -180/180 to a 0-360 grid.
shiftgrid_360_to_180(lon, data)	This function shifts ND data from a 0-360 to a -180/180 grid.
sol2ls(jld[, cumulative])	Return the solar longitude (Ls) as a function of the sol number.
sol_hhmmss(time_sol[, lon_180])	Given the time in days, return return local true solar time at a
spherical_curl(U, V, lon_deg, lat_deg[, R, spacing])	Compute the vertical component of the relative vorticity using
spherical_div(U, V, lon_deg, lat_deg[, R, spacing])	Compute the divergence of the wind fields using finite difference:
swinbank(plev, psfc[, ptrans])	Compute ak and bk values with a transition based on Swinbank
time_shift_calc(var_in, lon, tod[, target_times])	Conversion to uniform local time.
transition(pfull[, p_sigma, p_press])	Return the transition factor to construct ak and bk
vinterp(varIN, Lfull, Llev[, type_int, reverse_input, ...])	Vertical linear or logarithmic interpolation for pressure or altitude.
vw_from_MSF(msf, lat, lev[, ztype, norm, psfc, H])	Return the V and W components of the circulation from the mass
zonal_detrend(VAR)	Substract the zonal average mean value from a field.

MGStau_ls_lat(ls, lat)

Return the max altitude for the dust from “MGS scenario” from Montmessin et al. (2004), Origin and role of water ice clouds in the Martian water cycle as inferred from a general circulation model

Parameters:

ls (array) – solar longitude [°]

:param lat : latitude [°] :type lat: array :return: top altitude for the dust [km]

MGSzmax_ls_lat(ls, lat)

Return the max altitude for the dust from “MGS scenario” from Montmessin et al. (2004), Origin and role of water ice clouds in the Martian water cycle as inferred from a general circulation model

Parameters:

ls (array) – solar longitude [°]

:param lat : latitude [°] :type lat: array :return: top altitude for the dust [km]

UT_LTtxt(UT_sol, lon_180=0.0, roundmin=None)

Returns the time in HH:MM:SS at a certain longitude.

Parameters:

• time_sol (float) – the time in sols

• lon_180 (float) – the center longitude in -180/180 coordinates. Increments by 1hr every 15°

• roundmin (int) – round to the nearest X minute. Typical values are roundmin = 1, 15, 60

Note

If roundmin is requested, seconds are not shown

add_cyclic(data, lon)

Add a cyclic (overlapping) point to a 2D array. Useful for azimuth and orthographic projections.

Parameters:

• data (array) – variable of size [nlat, nlon]

• lon (array) – longitudes

Returns:

a 2D array of size [nlat, nlon+1] with last column identical to the 1st; and a 1D array of longitudes size [nlon+1] where the last element is lon[-1] + dlon

alt_KM(press, scale_height_KM=8.0, reference_press=610.0)

Gives the approximate altitude [km] for a given pressure

Parameters:

• press (1D array) – the pressure [Pa]

• scale_height_KM (float) – scale height [km] (default is 8 km, an isothermal at 155K)

• reference_press (float) – reference surface pressure [Pa] (default is 610 Pa)

Returns:

z_KM the equivalent altitude for that pressure [km]

Note

Scale height is H = rT/g

area_meridional_cells_deg(lat_c, dlon, dlat, normalize=False, R=3390000.0)

Return area of invidual cells for a meridional band of thickness dlon where S = int[R^2 dlon cos(lat) dlat] with sin(a)-sin(b) = 2 cos((a+b)/2)sin((a+b)/2) so S = 2 R^2 dlon 2cos(lat)sin(dlat/2):

    _________lat + dlat/2
    \    lat \               ^
    \lon +   \              | dlat
    \________\lat - dlat/2 v
lon - dlon/2   lon + dlon/2
        <------>
        dlon

Parameters:

• lat_c (float) – latitude of cell center [°]

• dlon (float) – cell angular width [°]

• dlat (float) – cell angular height [°]

• R (float) – planetary radius [m]

• normalize (bool) – if True, the sum of the output elements = 1

Returns:

S areas of the cells, same size as lat_c in [m2] or normalized by the total area

area_weights_deg(var_shape, lat_c, axis=-2)

Return weights for averaging the variable.

Expected dimensions are:

[lat] axis not needed [lat, lon] axis = -2 or axis = 0 [time, lat, lon] axis = -2 or axis = 1 [time, lev, lat, lon] axis = -2 or axis = 2 [time, time_of_day_24, lat, lon] axis = -2 or axis = 2 [time, time_of_day_24, lev, lat, lon] axis = -2 or axis = 3

Because dlat is computed as lat_c[1]-lat_c[0], lat_c may be truncated on either end (e.g., lat = [-20 ..., 0... 50]) but must be continuous.

Parameters:

• var_shape (tuple) – variable shape

• lat_c (float) – latitude of cell centers [°]

• axis (int) – position of the latitude axis for 2D and higher dimensional arrays. The default is the SECOND TO LAST dimension

Returns:

W weights for the variable ready for standard averaging as np.mean(var*W) [condensed form] or np.average(var, weights=W) [expanded form]

Note

Given a variable var:

var = [v1, v2, ...vn]

The regular average is

AVG = (v1 + v2 + ... vn) / N

and the weighted average is

AVG_W = (v1*w1 + v2*w2 + ... vn*wn) / (w1 + w2 + ...wn)

This function returns

W = [w1, w2, ... , wn]*N / (w1 + w2 + ...wn)

Therfore taking a regular average of (var*W) with np.mean(var*W) or np.average(var, weights=W)

returns the weighted average of the variable. Use np.average(var, weights=W, axis = X) to average over a specific axis.

areo_avg(VAR, areo, Ls_target, Ls_angle, symmetric=True)

Return a value average over a central solar longitude

EX:

``Ls_target = 90.``
``Ls_angle = 10.``

Nominally, the time average is done over solar longitudes 85 < Ls_target < 95 (10°).

If symmetric = True and the input data range is Ls = 88-100° then 88 < Ls_target < 92 (4°, symmetric)

If symmetric = False and the input data range is Ls = 88-100° then 88 < Ls_target < 95 (7°, assymetric)

Parameters:

• VAR (ND array) – a variable with time in the 1st dimension

• areo (1D array) – solar longitude of the input variable (0-720)

• Ls_target (float) – central solar longitude of interest

• Ls_angle (float) – requested window angle centered at Ls_target

• symmetric (bool (defaults to True)) – If True and the requested window is out of range, Ls_angle is reduced. If False, the time average is performed on the data available

Returns:

the variable averaged over solar longitudes Ls_target-Ls_angle/2 to Ls_target+Ls_angle/2

Note

The routine can bin data from muliples Mars years

axis_interp(var_IN, x, xi, axis, reverse_input=False, type_int='lin', modulo=None)

One dimensional linear/logarithmic interpolation along one axis.

Parameters:

• var_IN (ND array) – Variable on a REGULAR grid (e.g., [lev, lat, lon] or [time, lev, lat, lon])

• x (1D array) – Original position array (e.g., time)

• xi (1D array) – Target array to interpolate the array on

• axis (int) – Position of the interpolation axis (e.g., 0 for a temporal interpolation on [time, lev, lat, lon])

• reverse_input (bool) – Reverse input arrays (e.g., if ``zfull(0)``= 120 km, ``zfull(N)``= 0 km, which is typical)

• type_int (str) – “log” for logarithmic (typically pressure), “lin” for linear

• modulo (float) – For “lin” interpolation only, use cyclic input (e.g., when using modulo = 24 for time of day, 23.5 and 00 am are considered 30 min apart, not 23.5 hr apart)

Returns:

VAR_OUT interpolated data on the requested axis

Note

This routine is similar but simpler than the vertical interpolation vinterp() as the interpolation axis is assumed to be fully defined by a 1D array such as time, pstd or zstd rather than 3D arrays like pfull or zfull. For lon/lat interpolation, consider using interp_KDTree().

Calculation:

X_OUT = Xn*A + (1-A)*Xn + 1
with ``A = log(xi/xn + 1) / log(xn/xn + 1)`` in "log" mode
or ``A = (xi-xn + 1)/(xn-xn + 1)`` in "lin" mode

azimuth2cart(LAT, LON, lat0, lon0=0)

Azimuthal equidistant projection. Converts from latitude-longitude to cartesian coordinates.

Parameters:

• LAT (1D or 2D array) – latitudes[°] size [nlat]

• LON (1D or 2D array) – longitudes [°] size [nlon]

• lat0 (float) – latitude coordinate of the pole

• lon0 (float) – longitude coordinate of the pole

Returns:

the cartesian coordinates for the latitudes and longitudes

broadcast(var_1D, shape_out, axis)

Broadcast a 1D array based on a variable’s dimensions

Parameters:

• var_1D (1D array) – variable (e.g., lat size = 36, or time size = 133)

• shape_out (list) – broadcasting shape (e.g., temp.shape = [133, lev, 36, lon])

Returns:

(ND array) broadcasted variables (e.g., size = [1,36,1,1] for lat or [133,1,1,1] for time)

cart_to_azimut_TR(u, v, mode='from')

Convert cartesian coordinates or wind vectors to radians using azimuth angle.

Parameters:

• x (1D array) – the cartesian coordinate

• y (1D array) – the cartesian coordinate

• mode (str) – “to” for the direction that the vector is pointing, “from” for the direction from which the vector is coming

Returns:

Theta [°] and R the polar coordinates

compute_uneven_sigma(num_levels, N_scale_heights, surf_res, exponent, zero_top)

Construct an initial array of sigma based on the number of levels and an exponent

Parameters:

• num_levels (float) – the number of levels

• N_scale_heights (float) – the number of scale heights to the top of the model (e.g., N_scale_heights = 12.5 ~102 km assuming an 8 km scale height)

• surf_res (float) – the resolution at the surface

• exponent (float) – an exponent to increase the thickness of the levels

• zero_top (bool) – if True, force the top pressure boundary (in N = 0) to 0 Pa

Returns:

an array of sigma layers

daily_to_average(varIN, dt_in, nday=5, trim=True)

Bin a variable from an atmos_daily file format to the atmos_average file format.

Parameters:

• varIN (ND array) – variable with time dimension first (e.g., ts[time, lat, lon])

• dt_in (float) – delta of time betwen timesteps in sols (e.g., dt_in = time[1] - time[0])

• nday (int) – bining period in sols, default is 5 sols

• trim (bool) – whether to discard any leftover data at the end of file before binning

Returns:

the variable bin over nday

Note

If varIN[time, lat, lon] from atmos_daily is [160, 48, 96] and has 4 timesteps per day (every 6 hours), then the resulting variable for nday = 5 is varOUT(160/(4*5), 48, 96) = varOUT(8, 48, 96)

Note

If the daily file has 668 sols, then there are 133 x 5 + 3 sols leftover. If trim = True, then the time is 133 and last 3 sols the are discarded. If trim = False, the time is 134 and last bin contains only 3 sols of data.

daily_to_diurn(varIN, time_in)

Bin a variable from an atmos_daily file into the atmos_diurn format.

Parameters:

• varIN (ND array) – variable with time dimension first (e.g., [time, lat, lon])

• time_in (ND array) – time array in sols. Only the first N elements are actually required if saving memory is important

Returns:

the variable binned in the atmos_diurn format ([time, time_of_day, lat, lon]) and the time of day array [hr]

Note

If varIN[time, lat, lon] from atmos_daily is [40, 48, 96] and has 4 timestep per day (every 6 hours), then the resulting variable is varOUT[10, 4, 48, 96] = [time, time_of_day, lat, lon] and tod = [0., 6., 12., 18.].

Note

Since the time dimension is first, the output variables may be passed to the daily_to_average() function for further binning.

dvar_dh(arr, h=None)

Differentiate an array A[dim1, dim2, dim3...] w.r.t h. The differentiated dimension must be the first dimension.

EX: Compute dT/dz where T[time, lev, lat, lon] is the temperature and Zkm is the array of level heights [km].

First, transpose T so the vertical dimension comes first: T[lev, time, lat, lon].

Then transpose back to get dTdz[time, lev, lat, lon]:

dTdz = dvar_dh(t.transpose([1, 0, 2, 3]),
               Zkm).transpose([1, 0, 2, 3])

If h is 1D, then h``and ``dim1 must have the same length

If h is 2D, 3D or 4D, then arr and h must have the same shape

Parameters:

• arr (ND array) – variable

• h (str) – the dimension (Z, P, lat, lon)

Returns:

d_arr: the array differentiated w.r.t h, e.g., d(array)/dh

expand_index(Nindex, VAR_shape_axis_FIRST, axis_list)

Repeat interpolation indices along an axis.

Parameters:

• Nindex (idx) – Interpolation indices, size is (n_axis, Nfull = [time, lat, lon])

• VAR_shape_axis_FIRST (tuple) – Shape for the variable to interpolate with interpolation axis first (e.g., [tod, time, lev, lat, lon])

• axis_list (int or list) – Position or list of positions for axis to insert (e.g., 2 for lev in [tod, time, lev, lat, lon], [2, 4] for lev and lon). The axis positions are those for the final shape (VAR_shape_axis_FIRST) and must be INCREASING

Returns:

LFULL a 2D array (size n_axis, NfFULL = [time, lev, lat, lon]) with the indices expanded along the lev dimension and flattened

Note

Example of application: Observational time of day may be the same at all vertical levels so the interpolation of a 5D variable [tod, time, lev, lat, lon] only requires the interpolation indices for [tod, time, lat, lon]. This routine expands the indices from [tod, time, lat, lon] to [tod, time, lev, lat, lon] with Nfull = [time, lev, lat, lon] for use in interpolation.

find_n(X_IN, X_OUT, reverse_input=False, modulo=None)

Maps the closest index from a 1D input array to a ND output array just below the input values.

Parameters:

• X_IN (float or 1D array) – Source level [Pa] or [m]

• X_OUT (array) – Desired pressure [Pa] or altitude [m] at layer midpoints. Level dimension is FIRST

• reverse_input (bool) – If input array is decreasing (e.g., if z(0) = 120 km, z(N) = 0 km, which is typical, or if data is p(0) = 1000 Pa, p(N) = 0 Pa, which is uncommon)

Returns:

The index for the level(s) where the pressure < plev

find_n0(Lfull_IN, Llev_OUT, reverse_input=False)

Return the index for the level(s) just below Llev_OUT. This assumes Lfull_IN is increasing in the array (e.g., p(0) = 0, p(N) = 1000 [Pa]).

Parameters:

• Lfull_IN (array) – Input pressure [Pa] or altitude [m] at layer midpoints. Level dimension is FIRST

• Llev_OUT (float or 1D array) – Desired level type for interpolation [Pa] or [m]

• reverse_input (bool) – Reverse array (e.g., if z(0) = 120 km, z(N) = 0km – which is typical – or if input data is p(0) = 1000Pa, p(N) = 0Pa)

Returns:

n index for the level(s) where the pressure is just below plev

Note

If Lfull_IN is a 1D array and Llev_OUT is a float then n is a float.

Note

If Lfull_IN is ND [lev, time, lat, lon] and Llev_OUT is a 1D array of size klev then n is an array of size [klev, Ndim] with Ndim = [time, lat, lon]

fms_Z_calc(psfc, ak, bk, T, topo=0.0, lev_type='full')

Returns the 3D altitude field [m] AGL (or above aeroid).

Parameters:

• psfc (array) – The surface pressure [Pa] or array of surface pressures (1D, 2D, or 3D)

• ak (array) – 1st vertical coordinate parameter

• bk (array) – 2nd vertical coordinate parameter

• T (1D array or ND array) – The air temperature profile. 1D array (for a single grid point), ND array with VERTICAL AXIS FIRST

• topo (array) – The surface elevation. Same dimension as psfc. If None is provided, AGL is returned

• lev_type (str) – “full” (layer midpoint) or “half” (layer interfaces). Defaults to “full”

Returns:

The layer altitude at the full level Z_f(:, :, Nk-1) or half-level Z_h(:, :, Nk) [m]. Z_f and Z_h are AGL if topo = None. Z_f and Z_h are above aeroid if topography is not None.

Calculation:

--- 0 --- TOP        ========  z_half
--- 1 ---
                    --------  z_full

                    ========  z_half
---Nk-1---          --------  z_full
--- Nk --- SFC      ========  z_half
                    / / / / /

Note

Expands to the time dimension using:

topo = np.repeat(zsurf[np.newaxis, :], ps.shape[0], axis = 0)

Calculation is derived from

./atmos_cubed_sphere_mars/Mars_phys.F90:

# (dp/dz = -rho g) => (dz = dp/(-rho g)) and
# (rho = p/(r T)) => (dz = rT/g * (-dp/p))

# Define log-pressure (``u``) as:
u = ln(p)

# Then:
du = {du/dp}*dp = {1/p)*dp} = dp/p

# Finally, ``dz`` for the half-layers:
(dz = rT/g * -(du)) => (dz = rT/g * (+dp/p))
# with ``N`` layers defined from top to bottom.

Z_half calculation:

# Hydrostatic relation within the layer > (P(k+1)/P(k) =
# exp(-DZ(k)/H))

# layer thickness:
DZ(k) = rT/g * -(du)

# previous layer altitude + thickness of layer:
Z_h k) = Z_h(k+1)  +DZ_h(h)

Z_full calculation:

# previous altitude + half the thickness of previous layer and
# half of current layer
Z_f(k) = Z_f(k+1) + (0.5 DZ(k) + 0.5 DZ(k+1))

# Add ``+0.5 DZ(k)-0.5 DZ(k)=0`` and re-organiz the equation
Z_f(k) = Z_f(k+1) + DZ(k) + 0.5 (DZ(k+1) - DZ(k))
Z_f(k) = Z_h(k+1) + 0.5 (DZ(k+1) - DZ(k))

The specific heat ratio: γ = cp/cv (cv = cp-R) => γ = cp/(cp-R) Also (γ-1)/γ=R/cp

The dry adiabatic lapse rate: Γ = g/cp => Γ = (gγ)/R

The isentropic relation: T2 = T1(p2/p1)**(R/cp)

Therefore:

line 1) =====Thalf=====zhalf[k]          line 2)                                   line 3)                                    line 4) -----Tfull-----zfull[k]     \ T(z)= To-Γ (z-zo)
line 5)                                      line 6)                                       line 7) =====Thalf=====zhalf[k+1]

Line 1: T_half[k+1]/Tfull[k] = (p_half[k+1]/p_full[k])**(R/Cp)

Line 4: From the lapse rate, assume T decreases linearly within the layer so T_half[k+1] = T_full[k] + Γ(Z_full[k]-Z_half[k+1]) and (Tfull < Thalf and Γ > 0)

Line 7: Z_full[k] = Z_half[k] + (T_half[k+1]-T_full[k])/Γ Pulling out Tfull from above equation and using Γ = (gγ)/R:

Z_full[k] = (Z_half[k+1] + (R Tfull[k]) / (gγ)(T_half[k+1]
/ T_full[k] - 1))

Using the isentropic relation above:

Z_full = (Z_half[k+1] + (R Tfull[k]) / (gγ)(p_half[k+1]
/ p_full[k])**(R/Cp)-1))

fms_press_calc(psfc, ak, bk, lev_type='full')

Returns the 3D pressure field from the surface pressure and the ak/bk coefficients.

Parameters:

• psfc (array) – the surface pressure [Pa] or an array of surface pressures (1D, 2D, or 3D if time dimension)

• ak (array) – 1st vertical coordinate parameter

• bk (array:) – 2nd vertical coordinate parameter

• lev_type (str) – “full” (layer midpoints) or “half” (layer interfaces). Defaults to “full.”

Returns:

the 3D pressure field at the full levels PRESS_f(Nk-1:,:,:) or half-levels PRESS_h(Nk,:,:,) [Pa]

Calculation:

--- 0 --- TOP        ========  p_half
--- 1 ---
                     --------  p_full

                     ========  p_half
---Nk-1---           --------  p_full
--- Nk --- SFC       ========  p_half
                    / / / / /

Note

Some literature uses pk (pressure) instead of ak with p3d = ps * bk + P_ref * ak instead of p3d = ps * bk + ak

frontogenesis(U, V, theta, lon_deg, lat_deg, R=3400 * 1000.0, spacing='varying')

Compute the frontogenesis (local change in potential temperature gradient near a front) following Richter et al. 2010: Toward a Physically Based Gravity Wave Source Parameterization in a General Circulation Model, JAS 67.

We have Fn = 1/2 D(Del Theta)^2/Dt [K/m/s]

Parameters:

• U (array) – wind field with lat SECOND TO LAST and lon as last dimensions (e.g., [lat, lon] or [time, lev, lat, lon] etc.)

• V (array) – wind field with lat SECOND TO LAST and lon as last dimensions (e.g., [lat, lon] or [time, lev, lat, lon] etc.)

• theta (array) – potential temperature [K]

• lon_deg (1D array) – longitude [°] (2D if irregularly-spaced)

• lat_deg (1D array) – latitude [°] (2D if irregularly-spaced)

• R (float) – planetary radius [m]

• spacing (str (defaults to "varying")) – when lon and lat are 1D arrays, using spacing = “varying” differentiates latitude and longitude. When spacing = “regular”, dx = lon[1]-lon[0], `` dy=lat[1]-lat[0]`` and the numpy.gradient() method are used

Returns:

the frontogenesis field [m-1]

gauss_profile(x, alpha, x0=0.0)

Return Gaussian line shape at x. This can be used to generate a bell-shaped mountain.

get_trend_2D(VAR, LON, LAT, type_trend='wmean')

Extract spatial trends from the data. The output can be directly subtracted from the original field.

Parameters:

• VAR (ND array) – Variable for decomposition. lat is SECOND TO LAST and lon is LAST (e.g., [time, lat, lon] or [time, lev, lat, lon])

• LON (2D array) – longitude coordinates

• LAT (2D array) – latitude coordinates

• type_trend (str) – type of averaging to perform: “mean” - use a constant average over all lat/lon “wmean” - use a area-weighted average over all lat/lon “zonal” - detrend over the zonal axis only “2D” - use a 2D planar regression (not area-weighted)

Returns:

the trend, same size as VAR

interp_KDTree(var_IN, lat_IN, lon_IN, lat_OUT, lon_OUT, N_nearest=10)

Inverse distance-weighted interpolation using nearest neighboor for ND variables. Alex Kling, May 2021

Parameters:

• var_IN (ND array) – ND variable to regrid (e.g., [lev, lat, lon], [time, lev, lat, lon] with [lat, lon] dimensions LAST [°])

• lat_IN (1D or 2D array) – latitude [°] (LAT[y, x] array for irregular grids)

• lon_IN (1D or 2D array) – latitude [°] (LAT[y, x] array for irregular grids)

• lat_OUT (1D or 2D array) – latitude [°] for the TARGET grid structure (or LAT1[y,x] for irregular grids)

• lon_OUT (1D or 2D array) – longitude [°] for the TARGET grid structure (or LON1[y,x] for irregular grids)

• N_nearest (int) – number of nearest neighbours for the search

Returns:

VAR_OUT interpolated data on the target grid

Note

This implementation is much FASTER than griddata and it supports unstructured grids like an MGCM tile. The nearest neighbour interpolation is only done on the lon/lat axis (not level). Although this interpolation works well on the 3D field [x, y, z], this is typically not what is expected. In a 4°x4° run, the closest points in all directions (N, E, S, W) on the target grid are 100’s of km away while the closest points in the vertical are a few 10’s -100’s meter in the PBL. This would result in excessive weighting in the vertical.

layers_mid_point_to_boundary(pfull, sfc_val)

A general description for the layer boundaries is:

p_half = ps*bk + pk

This routine converts the coordinate of the layer MIDPOINTS, p_full or bk, into the coordinate of the layer BOUNDARIES p_half. The surface value must be provided.

Parameters:

• p_full (1D array) – Pressure/sigma values for the layer MIDPOINTS, INCREASING with N (e.g., [0.01 -> 720] or [0.001 -> 1])

• sfc_val (float) – The surface value for the lowest layer’s boundary p_half[N] (e.g., sfc_val = 720 Pa or sfc_val = 1 in sigma coordinates)

Returns:

p_half the pressure at the layer boundaries (size = N+1)

Structure:

--- 0 --- TOP   ========  p_half
--- 1 ---
                --------  p_full

                ========  p_half
---Nk-1---      --------  p_full
--- Nk --- SFC  ========  p_half
                / / / / /

We have:

pfull[N] = ((phalf[N]-phalf[N-1]) / np.log(phalf[N]/phalf[N-1]))
=> phalf[N-1] - pfull[N] log(phalf[N-1])
= phalf[N] - pfull[N] log(phalf[N])

We want to solve for phalf[N-1] = X:

v                v                             v
X      - pfull[N]       log(X)   =             B

=> X= -pfull[N] W{-exp(-B/pfull[N])/pfull[N]}

with B = phalf[N] - pfull[N] log(phalf[N]) (known at N) and

W is the product-log (Lambert) function.

This was tested on an L30 simulation: The values of phalf are reconstructed from pfull with a max error of:

100*(phalf - phalf_reconstruct)/phalf < 0.4% at the top.

lin_interp(X_in, X_ref, Y_ref)

Simple linear interpolation with no dependance on scipy

Parameters:

X_in (float or array) – input values

:param X_ref x values :type X_ref: array :param Y_ref y values :type Y_ref: array :return: y value linearly interpolated at X_in

lon180_to_360(lon)

Transform a float or an array from the -180/180 coordinate system to 0-360

Parameters:

lon (float, 1D array, or 2D array) – longitudes in the -180/180 coordinate system

Returns:

the equivalent longitudes in the 0-360 coordinate system

lon360_to_180(lon)

Transform a float or an array from the 0-360 coordinate system to

-180/180.

Parameters:

lon (float, 1D array, or 2D array) – longitudes in the 0-360 coordinate system

Returns:

the equivalent longitudes in the -180/180 coordinate system

ls2sol(Ls_in)

Ls to sol converter.

Parameters:

Ls_in (float or 1D array) – solar longitudes (0-360…720)

Returns:

the corresponding sol number

Note

This function simply uses a numerical solver on the sol2ls() function.

mass_stream(v_avg, lat, level, type='pstd', psfc=700, H=8000.0, factor=1e-08)

Compute the mass stream function:

                        P
                        ⌠
Ph i= (2 pi a) cos(lat)/g ⎮vz_tavg dp
                        ⌡
                        p_top

Parameters:

• v_avg (ND array) – zonal wind [m/s] with lev dimension FIRST and lat dimension SECOND (e.g., [pstd, lat], [pstd, lat, lon] or [pstd, lat, lon, time])

• lat (1D array) – latitudes [°]

• level (1D array) – interpolated layers [Pa] or [m]

• type (str) – interpolation type (pstd, zstd or zagl)

• psfc (float) – reference surface pressure [Pa]

• H (float) – reference scale height [m] when pressures are used

• factor (int) – normalize the mass stream function by a factor, use factor = 1 for [kg/s]

Returns:

MSF the meridional mass stream function (in factor * [kg/s])

Note

This routine allows the time and zonal averages to be computed before OR after the MSF calculation.

Note

The expressions for MSF use log(pressure) Z coordinates, which integrate better numerically.

With p = p_sfc exp(-Z/H) and Z = H log(p_sfc/p) then dp = -p_sfc/H exp(-Z/H) dZ and we have:

                                Z_top
                                ⌠
Phi = +(2pi a)cos(lat)psfc/(gH) ⎮v_rmv exp(-Z/H)dZ
                                ⌡
                                Z

With p = p_sfc exp(-Z/H)

The integral is calculated using trapezoidal rule:

    n
    ⌠
.g. ⌡ f(z)dz = (Zn-Zn-1){f(Zn) + f(Zn-1)}/2
  n-1

mollweide2cart(LAT, LON)

Mollweide projection. Converts from latitude-longitude to cartesian coordinates.

Parameters:

• LAT (1D or 2D array) – latitudes[°] size [nlat]

• LON (1D or 2D array) – longitudes [°] size [nlon]

• lat0 (float) – latitude coordinate of the pole

• lon0 (float) – longitude coordinate of the pole

Returns:

the cartesian coordinates for the latitudes and longitudes

ortho2cart(LAT, LON, lat0, lon0=0)

Orthographic projection. Converts from latitude-longitude to cartesian coordinates.

Parameters:

• LAT (1D or 2D array) – latitudes[°] size [nlat]

• LON (1D or 2D array) – longitudes [°] size [nlon]

• lat0 (float) – latitude coordinate of the pole

• lon0 (float) – longitude coordinate of the pole

Returns:

the cartesian coordinates for the latitudes and longitudes; and a mask (NaN array) that hides the back side of the planet

polar2XYZ(lon, lat, alt, Re=3400 * 10**3)

Spherical to cartesian coordinate transformation.

Parameters:

• lon (ND array) – Longitude in radians

• lat (ND array) – Latitude in radians

• alt (ND array) – Altitude [m]

• Re (float) – Planetary radius [m], defaults to 3400*10^3

Returns:

X, Y, Z in cartesian coordinates [m]

Note

This is a classic polar coordinate system with colatitude = pi/2 - lat where cos(colat) = sin(lat)

polar_warming(T, lat, outside_range=np.nan)

Return the polar warming, following McDunn et al. 2013: Characterization of middle-atmosphere polar warming at Mars, JGR Alex Kling

Parameters:

• T (ND array) – temperature with the lat dimension FIRST (transpose as needed)

• lat (1D array) – latitude array

• outside_range (float) – values to set the polar warming to when outside pf the range. Default = NaN but 0 may be desirable.

Returns:

The polar warming [K]

Note

polar_warming() concatenates the results from both hemispheres obtained from the nested function PW_half_hemisphere()

press_pa(alt_KM, scale_height_KM=8.0, reference_press=610.0)

Gives the approximate altitude [km] for a given pressure

Parameters:

• alt_KM (1D array) – the altitude [km]

• scale_height_KM (float) – scale height [km] (default is 8 km, an isothermal at 155K)

• reference_press (float) – reference surface pressure [Pa] (default is 610 Pa)

Returns:

press_pa the equivalent pressure at that altitude [Pa]

Note

Scale height is H = rT/g

press_to_alt_atmosphere_Mars(Pi)

Return the altitude [m] as a function of pressure from the analytical calculation above.

Parameters:

Pi (float or 1D array) – input pressure [Pa] (<= 610 Pa)

Returns:

the corresponding altitude [m] (float or 1D array)

ref_atmosphere_Mars_PTD(Zi)

Analytical atmospheric model for Martian pressure, temperature, and density. Alex Kling, June 2021

Parameters:

Zi (float or 1D array) – input altitude [m] (must be >= 0)

Returns:

tuple of corresponding pressure [Pa], temperature [K],

and density [kg/m3] floats or arrays

Note

This model was obtained by fitting globally and annually averaged reference temperature profiles derived from the Legacy GCM, MCS observations, and Mars Climate Database.

The temperature fit was constructed using quadratic temperature T(z) = T0 + gam(z-z0) + a*(z-z0)^2 over 4 segments (0>57 km, 57>110 km, 110>120 km and 120>300 km).

From the ground to 120 km, the pressure is obtained by integrating (analytically) the hydrostatic equation:

dp/dz=-g. p/(rT) with T(z) = T0 + gam(z-z0) + a*(z-z0)^2

Above ~120 km, P = P0 exp(-(z-z0)g/rT) is not a good approximation as the fluid is in molecula regime. For those altitudes, we provide a fit in the form of P = P0 exp(-az-bz^2) based on diurnal average of the MCD database at lat = 0, Ls = 150.

regression_2D(X, Y, VAR, order=1)

Linear and quadratic regression on the plane.

Parameters:

• X (2D array) – first coordinate

• Y (2D array) – second coordinate

• VAR (2D array) – variable of the same size as X

• order (int) – 1 (linear) or 2 (quadratic)

Note

When order = 1, the equation is: aX + bY + C = Z. When order = 2, the equation is: aX^2 + 2bX*Y + cY^2 + 2dX + 2eY + f = Z

For the linear case::, ax + by + c = z is re-written as A X = b with:

    |x0   y0   1|        |a      |z0
A = |x1   y1   1|    X = |b   b= |
    |      ...  |        |c      |...
    |xn   yn   1|                |zn

        [n,3]           [3]       [n]

The least-squares regression provides the solution that that minimizes ||b – A x||^2

robin2cart(LAT, LON)

Robinson projection. Converts from latitude-longitude to cartesian coordinates.

Parameters:

• LAT (1D or 2D array) – latitudes[°] size [nlat]

• LON (1D or 2D array) – longitudes [°] size [nlon]

• lat0 (float) – latitude coordinate of the pole

• lon0 (float) – longitude coordinate of the pole

Returns:

the cartesian coordinates for the latitudes and longitudes

second_hhmmss(seconds, lon_180=0.0)

Given the time [sec], return local true solar time at a certain longitude.

Parameters:

• seconds (float) – the time [sec]

• lon_180 (float) – the longitude in -180/180 coordinate

Returns:

the local time [float] or a tuple (hours, minutes, seconds)

sfc_area_deg(lon1, lon2, lat1, lat2, R=3390000.0)

Return the surface between two sets of latitudes/longitudes:

S = int[R^2 dlon cos(lat) dlat]     _____lat2
                                    \                                                 \____\lat1
                                     lon1    lon2

Parameters:

• lon1 (float) – longitude from set 1 [°]

• lon2 (float) – longitude from set 2 [°]

• lat1 (float) – latitude from set 1 [°]

• lat2 (float) – longitude from set 2 [°]

• R (int) – planetary radius [m]

qLon and Lat define the corners of the area not the grid cell center.

shiftgrid_180_to_360(lon, data)

This function shifts ND data from a -180/180 to a 0-360 grid.

Parameters:

• lon (1D array) – longitudes in the 0-360 coordinate system

• data (ND array) – variable with lon in the last dimension

Returns:

shifted data

shiftgrid_360_to_180(lon, data)

This function shifts ND data from a 0-360 to a -180/180 grid.

Parameters:

• lon (1D array) – longitudes in the 0-360 coordinate system

• data (ND array) – variable with lon in the last dimension

Returns:

shifted data

Note

Use np.ma.hstack instead of np.hstack to keep the masked array properties

sol2ls(jld, cumulative=False)

Return the solar longitude (Ls) as a function of the sol number. Sol=0 is the spring equinox.

Parameters:

• jld (float or 1D array) – sol number after perihelion

• cumulative (bool) – if True, result is cumulative (Ls=0-360, 360-720 etc..)

Returns:

the corresponding solar longitude

sol_hhmmss(time_sol, lon_180=0.0)

Given the time in days, return return local true solar time at a certain longitude.

Parameters:

• time_sol – the time in sols

• lon_180 (float) – the longitude in -180/180 coordinate

Returns:

the local time [float] or a tuple (hours, minutes, seconds)

spherical_curl(U, V, lon_deg, lat_deg, R=3400 * 1000.0, spacing='varying')

Compute the vertical component of the relative vorticity using finite difference:

curl = dv/dx -du/dy
     = 1/(r cos lat)[d(v)/dlon + d(u(cos lat)/dlat]

Parameters:

• U (array) – wind field with lat SECOND TO LAST and lon as last dimensions (e.g., [lat, lon] or [time, lev, lat, lon] etc.)

• V (array) – wind field with lat SECOND TO LAST and lon as last dimensions (e.g., [lat, lon] or [time, lev, lat, lon] etc.)

• lon_deg (1D array) – longitude [°] (2D if irregularly-spaced)

• lat_deg (1D array) – latitude [°] (2D if irregularly-spaced)

• R (float) – planetary radius [m]

• spacing (str (defaults to "varying")) – when lon and lat are 1D arrays, using spacing = “varying” differentiates latitude and longitude. When spacing = “regular”, dx = lon[1]-lon[0], `` dy=lat[1]-lat[0]`` and the numpy.gradient() method are used

Returns:

the vorticity of the wind field [m-1]

spherical_div(U, V, lon_deg, lat_deg, R=3400 * 1000.0, spacing='varying')

Compute the divergence of the wind fields using finite difference:

div = du/dx + dv/dy
-> = 1/(r cos lat)[d(u)/dlon + d(v cos lat)/dlat]

Parameters:

• U (array) – wind field with lat SECOND TO LAST and lon as last dimensions (e.g., [lat, lon] or [time, lev, lat, lon] etc.)

• V (array) – wind field with lat SECOND TO LAST and lon as last dimensions (e.g., [lat, lon] or [time, lev, lat, lon] etc.)

• lon_deg (1D array) – longitude [°] (2D if irregularly-spaced)

• lat_deg (1D array) – latitude [°] (2D if irregularly-spaced)

• R (float) – planetary radius [m]

• spacing (str (defaults to "varying")) – when lon and lat are 1D arrays, using spacing = “varying” differentiates latitude and longitude. When spacing = “regular”, dx = lon[1]-lon[0], `` dy=lat[1]-lat[0]`` and the numpy.gradient() method are used

Returns:

the horizonal divergence of the wind field [m-1]

swinbank(plev, psfc, ptrans=1.0)

Compute ak and bk values with a transition based on Swinbank

Parameters:

• plev (1D array) – the pressure levels [Pa]

• psfc (1D array) – the surface pressure [Pa]

• ptrans (1D array) – the transition pressure [Pa]

Returns:

the coefficients for the new layers

time_shift_calc(var_in, lon, tod, target_times=None)

Conversion to uniform local time.

Mars rotates approx. 14.6° lon per Mars-hour (360° ÷ 24.6 hr) Each 14.6° shift in lon represents a 1-hour shift in local time This code uses the more precise calculation: lon_shift = 24.0 * lon / 360.

Parameters:

• var_in (ND array) – variable to be shifted. Assume lon is the first dimension and tod is the last dimension

• lon (1D array) – longitude

• tod (1D array) – time_of_day index from the input file

• target_times (float (optional)) – local time(s) [hr] to shift to (e.g., "3. 15.")

Returns:

the array shifted to uniform local time

Note

If target_times is not specified, the file is interpolated on the same tod as the input

transition(pfull, p_sigma=0.1, p_press=0.05)

Return the transition factor to construct ak and bk

In the MGCM code, the full pressures are computed from:

              del(phalf)
pfull = -----------------------------
        log(phalf(k+1/2)/phalf(k-1/2))

Parameters:

• pfull (1D array) – the pressure [Pa]

• p_sigma (float) – the pressure level where the vertical grid starts transitioning from sigma to pressure

• p_press (float) – the pressure level above which the vertical grid is pure (constant) pressure

Returns:

the transition factor. = 1 for pure sigma, = 0 for pure pressure and =0-1 for the transition

vinterp(varIN, Lfull, Llev, type_int='log', reverse_input=False, masktop=True, index=None)

Vertical linear or logarithmic interpolation for pressure or altitude.

Parameters:

• varIN (ND array) – Variable to interpolate (VERTICAL AXIS FIRST)

• Lfull (array) – Pressure [Pa] or altitude [m] at full layers, same dimensions as varIN

• Llev (1D array) – Desired level for interpolation [Pa] or [m]. May be increasing or decreasing as the output levels are processed one at a time

• type_int (str) – “log” for logarithmic (typically pressure), “lin” for linear (typically altitude)

• reverse_input (bool) – Reverse input arrays. e.g., if zfull[0] = 120 km then zfull[N] = 0km (typical) or if input data is ``pfull[0]``=1000 Pa, ``pfull[N]``=0 Pa

• masktop (bool) – Set to NaN values if above the model top

• index (None or array) – Indices for the interpolation, already processed as [klev, Ndim]. Indices calculated if not provided

Returns:

varOUT variable interpolated on the Llev pressure or altitude levels

Note

This interpolation assumes pressure decreases with height:

--  0  -- TOP  [0 Pa]   : [120 km]| X_OUT = Xn*A + (1-A)*Xn + 1
--  1  --               :         |
                        :         |
--  n  -- pn   [30 Pa]  : [800 m] | Xn
                        :         |
--  k  -- Llev [100 Pa] : [500 m] | X_OUT
-- n+1 -- pn+1 [200 Pa] : [200 m] | Xn+1

-- SFC --
/ / / / / /

with A = log(Llev/pn + 1) / log(pn/pn + 1) in “log” mode or A = (zlev-zn + 1) / (zn-zn + 1) in “lin” mode

vw_from_MSF(msf, lat, lev, ztype='pstd', norm=True, psfc=700, H=8000.0)

Return the V and W components of the circulation from the mass stream function.

Parameters:

• msf (ND array) – the mass stream function with lev SECOND TO LAST and the lat dimension LAST (e.g., [lev, lat], [time, lev, lat], [time, lon, lev, lat])

• lat (1D array) – latitude [°]

• lev (1D array) – level [Pa] or [m] (pstd, zagl, zstd)

• ztype (str) – Use pstd for pressure so vertical differentation is done in log space.

• norm (bool) – if True, normalize lat and lev before differentiating to avoid having to rescale manually the vectors in quiver plots

• psfc (float) – surface pressure for pseudo-height when ztype = pstd

• H (float) – scale height for pseudo-height when ztype = pstd

Returns:

the meditional and altitude components of the mass stream function for plotting as a quiver or streamlines.

Note

The components are: [v]=  g/(2 pi cos(lat)) dphi/dz [w]= -g/(2 pi cos(lat)) dphi/dlat

zonal_detrend(VAR)

Substract the zonal average mean value from a field.

Parameters:

VAR (ND array) – variable with detrending dimension last

Returns:

detrented field (same size as input)

Note

RuntimeWarnings are expected if the slice contains only NaNs which is the case below the surface and above the model top in the interpolated files. This routine disables such warnings temporarily.