Hi

I would like to calculate the thermosteric sea level change for some on-going TerraFIRMA overshoot runs. I understand that the relevant CMIP variable is zostoga (monthly globally averaged) but that this was calculated from UKESM output post-run before submission to CMIP (i.e., is not standard UKESM output). Does anyone have code to do the calculation that they can share?

All the best

tony

Hi - found it in CDDS! `CDDS/mip_convert/build/lib/mip_convert/process/processors.py`

starting line 1107

```
def calc_zostoga(thetao, thkcello, areacello, zfullo_0, so_0, rho_0_mean, deptho_0_mean):
"""
Calculate the global mean thermosteric sea level change with
respect to a reference state.
Parameters
---------
thetao: :class:`iris.cube.Cube`
Potential temperature (degrees Celsius)
thkcello: :class:`iris.cube.Cube`
Cell thickness (m)
areacello: :class:`iris.cube.Cube`
Cell area (m2)
zfullo_0: :class:`iris.cube.Cube`
**Reference state** depth (m)
so_0: :class:`iris.cube.Cube`
**Reference state** practical salinity (PSU)
rho_0_mean: :class:`iris.cube.Cube`
**Reference state** global mean in-situ density (kg/m3)
deptho_0_mean: :class:`iris.cube.Cube`
**Reference state** global mean water column depth (m),
computed as the total volume (`volo`) divided by the total surface area
(sum of `areacello`)
Returns
-------
zostoga: :class:`iris.cube.Cube`
Global mean thermosteric sea level change (m)
Returns
-------
zostoga: :class:`iris.cube.Cube`
Global mean thermosteric sea level change (m)
Notes
-----
Following eq. H27 of Griffies et al. (2016) [calc_zostoga_1]_:
| `zostoga = (volo_0 / areao) * (1 - RHO_t)`
| `RHO_t = RHO(thetao, so_0, p_0) / rho_0_mean`
where:
| `RHO` = Global mean in-situ density (kg/m3)
| `p_0` = **Reference state** sea water pressure at model levels (dbar),
approximated by the depth `zfullo_0`
The reference state is taken to be the first annual mean from the
piControl experiment for a particular model configuration. This is used to
calculate zostoga for all simulations based on this configuration.
References
----------
.. [calc_zostoga_1]
Griffies, S.M., Danabasoglu, G., Durack, P.J., Adcroft, A.J.,
Balaji, V., Boning, C.W., Chassignet, E.P., Curchitser,
E., Deshayes, J., Drange, H. and Fox-Kemper, B.,
2016. OMIP contribution to CMIP6: experimental and diagnostic
protocol for the physical component of the Ocean Model
Intercomparison Project. Geoscientific Model Development,
9(9), pp.3231-3296.
"""
rho_mean = iris.cube.CubeList()
so_0 = iris.util.squeeze(so_0)
zfullo_0 = iris.util.squeeze(zfullo_0)
for t_slice, z_slice in zip(thetao.slices_over('time'), thkcello.slices_over('time')):
# Check that thetao and thkcello have the same time coordinate
t_time = t_slice.coord('time')
z_time = z_slice.coord('time')
if t_time.cell(0) != z_time.cell(0):
message = 'Time coordinates of thetao and thkcello do not match:\n\tthetao: {}\n\tthkcello: {}'
raise ValueError(message.format(t_time, z_time))
# Change in mean in-situ density due to temperature (vs reference state)
rho_mean += [calc_rho_mean(t_slice, so_0, zfullo_0, areacello, z_slice)]
rho_mean = rho_mean.merge_cube()
rho_t = rho_mean.data / rho_0_mean.data
# zostoga is calculated by restating the change in mean in-situ density as a change in mean sea surface height
zostoga = deptho_0_mean.data * (1. - rho_t)
zostoga = rho_mean.copy(zostoga)
zostoga.units = Unit('m')
return zostoga
```

All the best,

George