src.models package
Submodules
src.models.atmos module
src.models.atmos.
This model was initially written in python by Dr. Naomi Henderson.
It includes the atmospheric Matsuno-Gill model, and the Seager1991 surface flux model.
It was refactored by Simon Thomas into a class structure / to be pylint compatible, and to take the cfg DictConfig struct as input.
It includes both the tropical Atmospheric model (Matsumo-Gill) with a single baroclinic mode, and the surface fluxes calculated based on the Ocean model.
An introduction to the Matsumo-Gill model:
https://www.atmos.washington.edu/~dargan/591/591_8.pdf
pytest src/test/test_atmos.py python3 src/models/atmos.py
Model Solution Method from Seager et al. 2019 Appendix:
The atmosphere equations are solved by Fourier transforming in longitude, forming an equation for v for each zonal wavenumber that is finite differenced, and the resulting tri-diagonal system is solved by matrix inversion, transforming back into longitude. Finally, u and Φ are derived by back-substitution. The ocean equations are solved using the ‘INC’ scheme, integrating the model forward, after spin-up with climatological conditions, forced by the time-varying ECMWF wind stress and, for the case with CO2 forcing, changing 𝑓′1 in the net surface longwave radiation calculation. Change over 1958–2017 is computed by a linear trend. The atmosphere model is solved forced by a Ts comprised of the climatological mean for 1958–2017 plus and minus half of the SST trend and the difference of the two simulations taken to derive the change. For the coupled model, the ocean model is first forced with the change in CO2 and climatological wind stress over 1958–2017. The resulting SST trend, plus the imposed heating change over land, are used to force the atmosphere model. The ocean model is forced again with both the changed wind stress and the CO2 increase to derive a new SST change over 1958–2017 that is then used to force the atmosphere model. This iterative coupling is repeated until equilibrium is reached, which takes just a few times. There is a unique solution for any given value of CO2. The model wind stress change is computed as 𝜌a𝑐D𝑊𝐮, where cD is a drag coefficient and 𝐮 is the vector surface wind change computed by the atmosphere model, which is added to the ECMWF climatological stresses. Since the atmosphere model dynamics are only applicable in the tropics, the computed wind stress anomaly is only applied to the ocean model between 20 S and 20 N, and is linearly tapered to zero at 25 S and 2 N.
Example
Old code to record the initial parameter settings:
# ------------- constants -----------------------
# begining TCAM
self.pr_max: float = 20.0 / 3600 / 24
# 20 / seconds in hour / hours in day.
self.relative_humidity: float = 0.80 # relative humidity uniformly 0.8
self.number_iterations: int = 50 # int
self.R_earth.to_value() = 6.37e6 # metres
self.sec_in_day = 86400 # seconds in day.
self.omega_2 = 2 * (2 * np.pi / self.sec_in_day) # 2 * rad per second
self.latent_heat_vap = 2.5e6 # latent heat # J kg-1
self.cp_air = 1000 # cp_air is the specific heat capacity of air.
# J kg-1 K-1
self.b_coeff = (
self.gravity * np.pi / (self.nbsq * self.theta_00 * self.height_tropopause)
)
self.eps_days = 0.75
# Over 1958–2017, the CO2 changed from ~300 to ~400 ppm,
# which would be about 0.75 W m−2
# eps days might be the efficiency of entrainment.
self.eps = 1.0 / (self.eps_days * self.sec_in_day) # 1/.75 d
self.eps_u = self.eps # 1/.75 d
self.eps_v = self.e_frac * self.eps # e_frac=1/2 in paper
# Newtonian cooling, K
self.newtonian_cooling_coeff_k1 = self.b_coeff / (self.k_days * self.sec_in_day)
self.eps_p = (np.pi / self.height_tropopause) ** 2 / (
self.nbsq * self.k_days * self.sec_in_day
)
self.beta = self.omega_2 / self.R_earth.to_value()
self.rho_air = 1.225 # kg m-3 - also called rho_00
self.c_e = 0.00125 # 1.25e-3 # cE is an exchange coefficient
self.emmisivity = 0.97 # problem here with second definintion.
self.stefan_boltzman_const = 5.67e-8
self.p_s = 1000 # pressure at the surface?
self.es_0 = 6.11
self.delta_temp = 1.0 # ΔT = 1 K
self.f2 = 0.05 # f2 = 0.05
# 'a_cloud_const' should decrease when deep convection happens above 28 degC
# python main.py --multirun oc.nummode=6,7,8
# a_cloud_const = Ts-temp_0_c;
# a_cloud_const[a_cloud_const>28] = 40;
# a_cloud_const[a_cloud_const<=28] = 80;
# a_cloud_const = 0.01*a_cloud_const
self.a_cloud_const = 0.6 # this isn't the option used in the paper.
# basic parameters
self.temp_0_c = 273.15 # zero degrees in kelvin
self.f1_bar = 0.39 # f1 = 0.39
# f'1 is the anomaly in f1—a parameter that can be adjusted
# to control the variation in surface longwave radiation due
# to a_cloud_const change in CO2
self.u_bar = 5.0 # average velocity?
self.temp_surface_bar: float = self.temp_0_c + 25 # 25C in Kelvin
self.c_bar = 0.6 # C is the cloud cover. perhaps C_bar is the average.
- class src.models.atmos.Atmos(cfg, setup)
Bases:
objectAtmos class.
- f_cor(y_axis)
Corriolis force coeff. Makes beta plane approximation.
omega_2 = 2 * (2 * np.pi / sec_in_day) # 2 * rad per second
- Parameters
y (np.ndarray) – latitude
- Returns
Corriolis force coeff.
- Return type
np.ndarray
- f_dqlh_dtemp(temperature, u_sp, rh_loc)
Flux dqlh_dtemp.
- Parameters
temperature (xr.DataArray) – temperature (in kelvin).
u_sp (Union[xr.DataArray, float]) – u speed.
rh_loc (xr.DataArray) – relative humidity.
- Returns
flux dqlh_dtemp.
- Return type
xr.DataArray
- f_dqlw_df(temperature, cloud_cover)
Flux dqlw_df.
- Parameters
temperature (xr.DataArray) – temperature in Kelvin.
cloud_cover (Union[xr.DataArray, float]) – constant.
- Returns
flux dqlw_df.
- Return type
xr.DataArray
- f_dqlw_dtemp(temperature, cloud_cover, f, rh_loc)
Flux dqlw_dtemp.
- f_dqs_dtemp(temperature)
Flux dqs dtemp.
change in saturation humdity of the suface by temperature.
- Parameters
temperature (xr.DataArray) – temp.
- Returns
flux q_s.
- Return type
xr.DataArray
- f_ebar(temperature, rh_loc)
Flux e_bar.
- Parameters
temperature (xr.DataArray) – temperature in kelvin.
rh_loc (xr.DataArray) – Relative humidity (dimemsionless)
- Returns
evaporation.
- Return type
xr.DataArray
- f_es(temperature)
Flux es.
- Parameters
temperature (xr.DataArray) – temp.
- Returns
Flux es.
- Return type
xr.DataArray
- f_evap(mask, q_a, wnsp)
evaporation flux.
- Parameters
mask (np.ndarray) – land mask
q_a (np.ndarray) – surface air humidity
wnsp (np.ndarray) – surface windspeed in m/s
- Returns
Evap in kg/m^2/s.
- Return type
np.ndarray
- f_mc(q_a, u, v)
moisture convergence flux.
To calculate surface heat fluxes and atmospheric moisture convergence, relative humidity is assumed to be spatially uniform in our standard model. (N.B., v is on y_axis_v points, u,q are on y_axis_u points)
- Parameters
q_a (np.ndarray) – surface air humidity
u (np.ndarray) – low level winds in m/s
v (np.ndarray) – low level winds in m/s
- Returns
Moisture convergence in kg/m^2/s.
- Return type
np.ndarray
- f_qa(t_s, s_p)
Flux qa.
- Parameters
t_s (np.ndarray) – sst in Kelvin.
s_p (np.ndarray) – surface pressure in mb.
- Returns
q_a, surface specific humidity.
- Return type
np.ndarray
- f_qa2(temp_surface)
flux qa2.
- Parameters
temp_surface (np.ndarray) – sst in Kelvin.
- Returns
q_s, surface specific humidity.
- Return type
np.ndarray
- f_qlh(temperature, u_sp, rh_loc)
Heat flux from latent heat.
It is assumed that the surface heat flux anomaly is dominated by longwave and latent heat fluxes and that the solar radiation does not change and the sensible heat flux anomaly is small.
- Parameters
temperature (xr.DataArray) – the temperature dataarray.
u_sp (Union[xr.DataArray, float]) – the wind speed.
rh_loc (xr.DataArray) – relative humidity
- Returns
flux qlh.
- Return type
xr.DataArray
- f_qlw(temp, cloud_cover, f, rh_loc)
Full long wave heat flux.
- Parameters
temp (xr.DataArray) – temperature in kelvin.
cloud_cover (xr.DataArray) – cloud cover, float.
f (float) –
rh_loc (xr.DataArray) – relative humidity.
- Returns
Full long wave heat flux.
- Return type
xr.DataArray
- f_qlw1(temperature, cloud_cover, f, rh_loc)
The first term of the long wave flux equation (14).
Qlw1 = epsilon sigma T^4 f’ (1 - a_cloud_const C^2)
- Parameters
- Returns
The first term of the long wave flux equation.
- Return type
xr.DataArray
- f_qlw2(temperature)
Second term in QLW equation (14).
- Parameters
temperature (xr.DataArray) – temperature in Kelvin.
- Returns
heat from long wave heat flux term 2.
- Return type
xr.DataArray
- f_qs(temperature)
Flux q_s.
q_s(Ts) is the saturation-specific humidity at the SST rqs(Ts) = q_s(Ts)
- Parameters
temperature (xr.DataArray) – temp.
- Returns
Flux q_s.
- Return type
xr.DataArray
- f_temp_a(temperature)
Temperature anomaly.
Delta temp = 1 in paper.
- Parameters
temperature (xr.DataArray) – temperature in kelvin.
- Returns
temperature anomaly.
- Return type
xr.DataArray
- get_cloud_const(temperature)
Function to produce variable cloud constant using temperature.
Models an approximation to the effect of deep convection.
- Parameters
temperature (xr.DataArray) – temperature in Kelvin.
- Returns
- The a constant. Normally 0.4 or 0.8. A standard value of
0.6 is applied if you disable deep convection.
- Return type
xr.DataArray
- get_dclim()
Opens the files, and applies functions to get surface fluxes.
Get the surface fluxes, qd_df, dq_dT among other things.
Example
Logic behind calculating the flux parameters:
dq_dt = self.f_dqlh_dtemp(t_sb, u_b, rh_b) + self.f_dqlw_dtemp(t_sb, c_b, self.atm.f1_bar, rh) dq_df = self.f_dqlw_df(t_sb, c_b) dq_dt = (${atm.rho_air} * ${atm.c_e} * ${atm.latent_heat_vap}* u_sp * (self.atm.e_factor * self.atm.es_0 * np.exp(17.67 * (temperature - self.atm.temp_0_c) / (temperature - self.atm.temp_0_c + 243.5) ) / self.atm.p_s * (17.67 * 243.5) / (temperature - self.atm.temp_0_c + 243.5) ** 2 ) * (1 - rh) + ${atm.emmisivity} * ${atm.stefan_boltzman_const} * ( (1 - a_cloud_const * cloud_cover ** 2) * temperature ** 3 * ( 4 * f - self.atm.f2 * np.sqrt(e_bar) * (4 + temperature * dqs_dtemp / 2 / q_s) ) + 12 * temperature ** 2 * self.atm.delta_temp ) ) dq_df = ( ${atm.emmisivity} * ${atm.stefan_boltzman_const} * (1 - a_cloud_const * cloud_cover ** 2) * temperature ** 4 )- Return type
Tuple[Dataset,DataArray,DataArray,DataArray,DataArray,DataArray,DataArray,DataArray,DataArray]- Returns
- Tuple[
xr.Dataset, xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray,
]: dclim, u_b, alh, alw, blw, dtemp_se, rh, c_b, t_sb.
- load_clim60()
Load the inputs to get_dclim before processing.
This is ok for the first iteration, but will need to be changed.
- Returns
An mfdataset with “ts”, “clt”, “sfcWind”, “rh”.
- Return type
xr.Dataset
- make_figure(cmap='viridis', lat='latitude', lon='longitude')
Make figure.
- output_trends()
Output trends ds.
Runs the Matsuno-Gill model with the trends in preipitation and sea surface temperature half added and half subtracted to work out the change in the other variables.
𝜀𝑢 . 𝑢 − 𝑓 . 𝑣 + 𝜙 . 𝑥 = 0 (1)
𝜀𝑣 . 𝑣 + 𝑓 . 𝑢 + 𝜙 . 𝑦 = 0 (2)
𝜀𝜙 . 𝜙 + 𝑢 . 𝑥 + 𝑣 . 𝑦 = −𝑄1 (3)
- Return type
- s91_solver(q1)
S91 solver from TCAM.py.
Presumably name refers to Richard Seager 1991.
‘s91_solver’ 0.00564 s
The atmosphere equations are solved by Fourier transforming in longitude, forming an equation for v for each zonal wavenumber that is finite differenced, and the resulting tri-diagonal system is solved by matrix inversion, transforming back into longitude.
Usef fft, ifft.
g . pi . N ^ 2
- q1 = —————- . (k . theta_s . Q_c)
theta_00 . z_t
- Parameters
q1 (np.ndarray) – modified heating that drives winds.
- Returns
u, v, phi
- Return type
Tuple[np.ndarray, np.ndarray, np.ndarray]
- smooth121(da, sdims, number_smooths=1, perdims=[])
Applies [0.25, 0.5, 0.25] stencil in sdims, one at a time.
- Parameters
da (xr.DataArray) – xarray.DataArray - e.g., ds.var
sdims (list) – list of dimensions over which to smooth - e.g., [‘lat’,’lon’]
number_smooths (int, optional) – integer number of smooths to apply - e.g., 1. Defaults to 1.
perdims (list, optional) – list of dimension to be treated as period boundaries - e.g., [‘lon’]. Defaults to [].
- Returns
smoothed output.
- Return type
xr.DataArray
- tdma_solver(ny_loc, a_loc, b, c, d)
tdma solver.
‘tdma_solver’ 0.00243 s
Tri Diagonal Matrix Algorithm (a.k.a. Thomas algorithm) solver
E.g. https://gist.github.com/cbellei/8ab3ab8551b8dfc8b081c518ccd9ada9
- Parameters
ny_loc (int) – local version of ny_loc.
a_loc (np.ndarray) – [description]
b (np.ndarray) – [description]
c (np.ndarray) – [description]
d (np.ndarray) – [description]
- Returns
xc
- Return type
np.ndarray
src.models.coupling module
Coupling between ocean and atmospheric models.
Example
Import statement usage:
from src.models.coupling import Coupling
- class src.models.coupling.Coupling(cfg, setup)
Bases:
objectCoupled model part.
Change over 1958–2017 is computed by a linear trend. The atmosphere model is solved forced by a Ts comprised of the climatological mean for 1958–2017 plus and minus half of the SST trend and the difference of the two simulations taken to derive the change. For the coupled model, the ocean model is first forced with the change in CO2 and climatological wind stress over 1958–2017. The resulting SST trend, plus the imposed heating change over land, are used to force the atmosphere model. The ocean model is forced again with both the changed wind stress and the CO2 increase to derive a new SST change over 1958–2017 that is then used to force the atmosphere model. This iterative coupling is repeated until equilibrium is reached, which takes just a few times. There is a unique solution for any given value of CO2. The model wind stress change is computed as ρcWu*cD, where cD is a drag coefficient and u is the vector surface wind change computed by the atmosphere model, which is added to the ECMWF climatological stresses. Since the atmosphere model dynamics are only applicable in the tropics, the computed wind stress anomaly is only applied to the ocean model between 20S and 20N, and is linearly tapered to zero at 25S and 25N.
- f_stress(wind_speed_mean, u_wind, v_wind)
Wind stress flux.
\begin{equation} \vec{\tau}= \rho c_{D} W \vec{u} \end{equation}- Parameters
wind_speed_mean (Union[xr.DataArray, float]) – W, the climatological annual mean wind speed, which is taken from ECMWF reanalysis for our standard model and from the CMIP5 multimodel mean when examining causes of bias in the CMIP5 model.
u_wind (xr.DataArray) – u_u, the zonal wind velocity.
v_wind (xr.DataArray) – u_v, the meridional wind velocity.
- Returns
- tau_x (zonal wind stress),
tau_y (meridional wind stress).
- Return type
Tuple[xr.DataArray, xr.DataArray]
- get_tau_anom(wind, u_vel, v_vel)
Return the tau anomaly with a clipping.
- Parameters
wind (xr.DataArray) – wind speed field.
u_vel (xr.DataArray) – u wind velocity. (X, Yu).
v_vel (xr.DataArray) – v wind velocity. (X, Yv).
- Returns
tau_u, tau_v
- Return type
Tuple[xr.DataArray, xr.DataArray]
- log(it)
Log the important information about the run.
- replace_stress(it)
Replace the stress files.
Currently just resaves the clim files with a diff name.
- TODO: systematically explore other options:
linear change between taubeg and tauend
take half trend away or not.
turn off land precipitation.
make deep convection universal.
change in thermocline opposite of that expected - is this a sign error?
- replace_surface_temp(it)
Replace sst for forcing atmosphere model.
TODO: replace masking with actual mask.
- tau_anom_ds()
Wind stress anomaly. # TODO: Could simplify this function, and the following function.
- Returns
dataset with different different tau fields.
- Return type
xr.Dataset
src.models.model_setup module
Set up the model, copy the files, get the names.
- class src.models.model_setup.ModelSetup(direc, cfg, make_move=True)
Bases:
objectInitialise, store, and setup directories for models.
- stress_file()
# these are the default start files. # most of these are changed during the run # TODO remove these file names, move control elsewhere. wind_clim_file: tau-ECMWF-clim wind_file: tau-ECMWF dq_dtemp_file: dQdT-sample.nc dq_df_file: dQdf-sample.nc sst_file: sst-ECMWF-clim.nc mask_file: om_mask.nc
- Return type
src.models.ocean module
Ocean model.
- class src.models.ocean.Ocean(cfg, setup)
Bases:
objectOcean model component.
- copy_old_io(it)
Move old io files to their new name after that step has completed.
- edit_inputs(it)
Edit the input files for each iteration.
- run(command)
Runs a line of bash in the ocean/RUN directory and times how long it takes.
- src.models.ocean.replace_item(init, fin, loc_string_list)
Replace items in a list of strings.
A helper function for the Ocean class, as it allows the inputs to fortran/C ocean model to be rewritten.
src.models.poly module
To apply polynomial fits, and propogate error.
- src.models.poly.fit(x_npa, y_npa, reg_type='lin')
Fit a polynomial curve, with an estimate of the uncertainty.
- Parameters
- Returns
- Paramaters with uncertainty,
function to put data into.
- Return type
Tuple[unp.uarray, Callable]
- src.models.poly.plot(x_values, y_values, reg_type='lin', x_label='x label', y_label='y label', fig_path=None, ax_format='both')
Plot the polynomial.
- Parameters
x_values (Sequence[Union[float, int]]) – The x values to fit.
y_values (Sequence[Union[float, int]]) – The y values to fit.
reg_type (str, optional) – Which regression to do. Defaults to “lin”.
x_label (str, optional) – X label for plot. e.g. r”$Delta bar{T}_s$ over tropical pacific (pac) region [$Delta$ K]”
y_label (str) – Y labelsfor plot. e.g. r”$Delta bar{T}_s$ over nino3.4 region [$Delta$ K]”
fig_path (Optional[str], optional) – Path to stor the figure in. Defaults to None.
ax_format (Literal["both", "x", "y"], optional) – which axes to format in scientific notation. Defaults to “both”.
- Returns
- Paramaters with uncertainty,
function to put data into.
- Return type
Tuple[unp.uarray, Callable]
Example
Using plot:
import wandb_summarizer.download from src.models.poly import plot run_info = wandb_summarizer.download.get_results("sdat2/seager19") f_df = pd.DataFrame(run_info)[3:13] f_df = f_df.drop(labels=[11], axis=0) nino_3_list = f_df["end_trend_nino3"].tolist() nino_4_list = f_df["end_trend_nino4"].tolist() plot( nino_3_list, nino_4_list, x_label=r"$\Delta \bar{T}_s$ over nino3 region [$\Delta$ K]", y_label=r"$\Delta \bar{T}_s$ over nino4 region [$\Delta$ K]", ax_format=None, reg_type="lin", )