Code for zostoga

1

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

2

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