""" The hybrid sigma-pressure coordinate: where a model level actually is. Datasets are built in memory, the way the rest of the unit tests do it, with one case run against the real start.nc to keep the arithmetic honest against a file LMDZ actually wrote. """ import numpy as np import pytest import xarray as xr from conftest import open_sample from dispnc import paths, vertical def a_hybrid_ds(nlev=6, ps=None, altitude_km=True, mids=True, edges=True): """ A minimal LMDZ-shaped file: hybrid coefficients, a surface pressure, and the pseudo-altitude they imply at 610 Pa with a 10 km scale height. """ # Sigma at the ground, pressure at the top, as every model does it b_edge = np.linspace(1.0, 0.0, nlev + 1) ** 2 a_edge = 610.0 * (np.linspace(1.0, 0.0, nlev + 1) - b_edge) * 0.5 a_mid, b_mid = 0.5 * (a_edge[:-1] + a_edge[1:]), 0.5 * (b_edge[:-1] + b_edge[1:]) reference = a_mid + b_mid * 610.0 altitude = -10000.0 * np.log(reference / 610.0) if ps is None: ps = np.full((3, 4), 610.0) data = {} if mids: data['aps'] = (('altitude',), a_mid, {'units': 'Pa'}) data['bps'] = (('altitude',), b_mid, {}) if edges: data['ap'] = (('interlayer',), a_edge, {'units': 'Pa'}) data['bp'] = (('interlayer',), b_edge, {}) data['ps'] = (('latitude', 'longitude'), np.asarray(ps, dtype=float), {'units': 'Pa'}) return xr.Dataset( data, coords={ 'altitude': ('altitude', altitude / (1000.0 if altitude_km else 1.0), {'units': 'km' if altitude_km else 'm'}), 'latitude': ('latitude', np.linspace(60, -60, np.shape(ps)[0])), 'longitude': ('longitude', np.linspace(-180, 90, np.shape(ps)[1])), }) class Axis: def __init__(self, values, units=None, dim='altitude'): self.values, self.units, self.dim = np.asarray(values), units, dim @property def size(self): return self.values.size class Plan: kind = 'globe' def __init__(self, ds): self.z = Axis(ds['altitude'].values, ds['altitude'].attrs.get('units')) self.y = Axis(ds['latitude'].values, dim='latitude') self.x = Axis(ds['longitude'].values, dim='longitude') # --- reading the coefficients ------------------------------------------------ def test_mid_layers_are_rebuilt_from_the_interfaces(): """ A file may carry either set. The mid-layer value LMDZ writes is the mean of the two interfaces around it, so one set gives the other. """ ds = a_hybrid_ds() both = vertical.mid_coefficients(ds, 6) from_edges = vertical.mid_coefficients(ds.drop_vars(['aps', 'bps']), 6) assert np.allclose(both[0], from_edges[0]) assert np.allclose(both[1], from_edges[1]) def test_interfaces_are_rebuilt_from_the_mid_layers(): """ Rebuilt interfaces have to hold the whole column: the top one is pinned to zero pressure, as it is in every file LMDZ writes, or the mass integral would quietly lose the layer above the highest level. """ ds = a_hybrid_ds() rebuilt = vertical.interface_coefficients(ds.drop_vars(['ap', 'bp']), 6) assert rebuilt is not None a_edge, b_edge = rebuilt assert a_edge.size == b_edge.size == 7 assert (a_edge[-1], b_edge[-1]) == (0.0, 0.0) assert b_edge[0] == pytest.approx(1.0, abs=0.05), "the bottom is the ground" def test_the_names_are_matched_whatever_their_case(): """ CMOR and CESM write `PS`, `hyai`/`hybi`; the difference between `ps` and `PS` is not worth falling back to pseudo-altitudes over. """ ds = a_hybrid_ds().rename({'ps': 'PS', 'aps': 'HYAM', 'bps': 'HYBM', 'ap': 'HYAI', 'bp': 'HYBI'}) assert vertical.surface_pressure_name(ds) == 'PS' assert vertical.mid_coefficients(ds, 6) is not None assert vertical.interface_coefficients(ds, 6) is not None def test_a_file_without_the_coordinate_gives_nothing(): """ No coefficients means no better answer than the vertical axis itself, and the caller has to be told so rather than handed a guess. """ ds = a_hybrid_ds().drop_vars(['ap', 'bp', 'aps', 'bps']) assert vertical.air_altitudes(ds, Plan(ds), ds['ps'].values) is None def test_a_surface_pressure_on_another_grid_is_refused(): """ A file can hold a surface pressure on a staggered or coarser grid. Broadcasting that into place would put the levels somewhere plausible and wrong, so it is refused and the vertical axis is used as it stands. """ ds = a_hybrid_ds() assert vertical.air_altitudes(ds, Plan(ds), np.full((5, 7), 610.0)) is None def test_the_level_count_has_to_match(): """ A field on the interlayers is not described by the mid-layer coefficients, and pairing them anyway would place every level one half-layer wrong. """ ds = a_hybrid_ds(nlev=6) plan = Plan(ds) plan.z = Axis(np.arange(7.0), 'km') assert vertical.air_altitudes(ds, plan, ds['ps'].values) is None # --- gravity ----------------------------------------------------------------- def test_gravity_comes_from_the_file_when_it_says_so(): """ LMDZ writes its constants into `controle`, so a run of the generic model on another planet is weighed with that planet's gravity. """ ds = a_hybrid_ds() assert vertical.gravity(ds) == pytest.approx(paths.planet_gravity) control = np.zeros(100) control[6] = 9.81 with_control = ds.assign(controle=('index', control)) assert vertical.gravity(with_control) == pytest.approx(9.81) def test_a_nonsense_gravity_is_ignored(): """ An empty or zeroed `controle` - a file written by something that is not LMDZ - must not divide the column by nothing. """ ds = a_hybrid_ds().assign(controle=('index', np.zeros(100))) assert vertical.gravity(ds) == pytest.approx(paths.planet_gravity) # --- the scale height -------------------------------------------------------- def test_the_scale_height_is_fitted_to_the_file_own_vertical(): """ `altitude` is -H*ln(p/p_ref) over a reference column, so H is recoverable from the file and never has to be assumed. It matters: the same tool is pointed at the generic model, whose atmosphere is not Mars'. """ ds = a_hybrid_ds() reference = ds['aps'].values + ds['bps'].values * 610.0 height, note = vertical.fit_scale_height(ds['altitude'].values, 'km', reference) assert height == pytest.approx(10000.0, rel=1e-6) assert note == '' def test_a_vertical_that_is_not_a_length_falls_back_on_the_constant(): """ Pressure or sigma levels give nothing to fit, and the caller is told which scale height it ended up with. """ height, note = vertical.fit_scale_height(np.array([1000.0, 500.0, 100.0]), 'Pa', np.array([1000.0, 500.0, 100.0])) assert height == pytest.approx(paths.scale_height) assert 'scale height' in note # --- the altitudes themselves ------------------------------------------------ def test_the_levels_follow_the_terrain_and_flatten_out_above_it(): """ The whole point of the hybrid coordinate: at the bottom B is one and the level *is* the ground, at the top A carries the pressure and the level is an isobar. Over a column with 25% less surface pressure - a mountain about a fifth of a scale height high - the lowest level still sits on the ground while the highest has to end up at the same absolute altitude as its neighbour, which means lower above the ground by exactly the height of the mountain. """ ps = np.full((3, 4), 610.0) ps[1, 2] = 610.0 * 0.75 # a summit, in surface pressure ds = a_hybrid_ds(nlev=20, ps=ps) air = vertical.air_altitudes(ds, Plan(ds), ps) assert air is not None assert np.all(np.diff(air.above_ground, axis=0) > 0), "levels must stack upwards" # The lowest level of twenty is at 0.95 ps, half a kilometre up: the model's # own first layer, not a shell floating over the terrain assert air.above_ground[0].max() < 600.0, "the lowest level is the ground" mountain = 10000.0 * np.log(1.0 / 0.75) # H * ln(ps0/ps), about 2.9 km over_summit = air.above_ground[-1, 1, 2] over_plain = air.above_ground[-1, 0, 0] # The top level gives back almost the whole of the mountain, so that the two # columns reach the same absolute altitude: the level is an isobar, and an # isobar does not care what is underneath it. Almost, not exactly - the # coefficient B is small at the top rather than zero, so a trace of the # terrain survives all the way up, which is what a hybrid coordinate is. assert over_plain - over_summit == pytest.approx(mountain, rel=0.1) absolute_summit = over_summit + mountain assert absolute_summit == pytest.approx(over_plain, rel=0.01) def test_the_thickness_of_the_whole_column_is_the_surface_pressure(): """ sum(|dp|)/g over the column is ps/g by construction, which is the check that the interfaces line up with the levels and nothing is counted twice. """ ps = np.linspace(500.0, 700.0, 12).reshape(3, 4) ds = a_hybrid_ds(ps=ps) air = vertical.air_altitudes(ds, Plan(ds), ps) assert np.allclose(air.thickness().sum(axis=0), ps / air.gravity) def test_a_column_with_no_surface_pressure_is_filled_rather_than_dropped(): """ A gap in `ps` would put a whole column of the shell nowhere in particular. The planetary mean is the least surprising stand-in, and the levels there still stack. """ ps = np.full((3, 4), 610.0) ps[0, 0] = np.nan ds = a_hybrid_ds(ps=ps) air = vertical.air_altitudes(ds, Plan(ds), ps) assert np.all(np.isfinite(air.above_ground)) assert air.above_ground[-1, 0, 0] == pytest.approx(air.above_ground[-1, 1, 1]) # --- against a file LMDZ really wrote ---------------------------------------- def test_start_nc_reproduces_its_own_altitude_axis(): """ Run on the real thing: the altitudes rebuilt from `aps`, `bps` and `ps` must come back to the file's own pseudo-altitude, which is what says the scale height, the coefficients and the reference pressure all agree. """ with open_sample('start.nc', 'ps', 'aps', 'bps', 'altitude') as (ds, _): plan = Plan(ds) air = vertical.air_altitudes(ds, plan, ds['ps'].values[0]) assert air is not None assert air.gravity == pytest.approx(3.72, abs=0.01), "Mars, from controle" published = ds['altitude'].values * 1000.0 profile = air.above_ground.reshape(published.size, -1).mean(axis=1) # Within a few per cent over the whole column: the mean is over every # surface pressure on the planet, rather than over the one reference # column the file's own axis was written from assert np.allclose(profile, published, rtol=0.05, atol=50.0) assert np.allclose(air.thickness().sum(axis=0), ds['ps'].values[0] / air.gravity, rtol=1e-6)