Verify diffusion model restart files and provide documentation to create from `.npy` files

[ma595]

Below is a very simple MWE to verify the metric results from spinup-evaluation. The outputs below agree with our spinup-evaluation results.

• depth_box == 3000 not 1500
• Average over volume not area
• Promote x (deg) and y (deg) to dimensions in restart file

TODO:

• Confirm the result for temperature_BWbox_metric
• Confirm for DWbox

import xarray as xr

path_to_restart_ref = "../spinup-data/Gorce-data/restart_files_1_4deg/reference_restart_file/DINO_12960000_restart.nc"

# path_to_restart_ref_david = "../spinup-data/Gorce-data/restart_files_1_4deg/restart_C2_David/DINO_11232000_restart.nc"

restart_ref = xr.load_dataset(path_to_restart_ref)
# restart_ref_david = xr.load_dataset(path_to_restart_ref_david)

mesh_mask = xr.load_dataset("../spinup-data/Gorce-data/Dinonline/restart0/mesh_mask.nc")


def temperature_BWbox_metric(thetao, file_mask, depth_box=1500):
    t_BW = thetao.where(1 - (thetao.nav_lev < depth_box) * (abs(thetao.y) < 30))

    # Computing Area Weights from Mask over Box
    e1t = file_mask.e1t.squeeze()
    e2t = file_mask.e2t.squeeze()
    e3t = file_mask.e3t_0.squeeze()
    tmask = file_mask.tmask.squeeze()
    area_BW = (
        e1t
        * e2t
        * e3t
        * tmask.where(1 - (thetao.nav_lev < depth_box) * (abs(thetao.y) < 30))
    )

    # Returning Average Temperature on Box
    return (t_BW * area_BW).sum(dim=["y", "x", "nav_lev"]) / area_BW.sum(
        dim=["y", "x", "nav_lev"]
    )


def temperature_DWbox_metric(thetao, file_mask, depth_box=1500):
    """
    Compute the average temperature in the deep water box defined as the region
    """

    e1t = file_mask.e1t.squeeze()
    e2t = file_mask.e2t.squeeze()
    e3t = file_mask.e3t_0.squeeze()
    tmask = file_mask.tmask.squeeze()
    condition = (
        (thetao.nav_lev < depth_box) * (thetao.nav_lev > 500) * (abs(thetao.y) < 30)
    )
    t_DW = thetao.where(condition)

    # Computing Area Weights from Mask over Box
    area_DW = e1t * e2t * e3t * tmask.where(condition)

    # Returning Average Temperature on Box
    return (t_DW * area_DW).sum(dim=["y", "x", "nav_lev"]) / area_DW.sum(
        dim=["y", "x", "nav_lev"]
    )

temperature = restart_ref.tb

print(
    "Bottom water temperature:",
    temperature_BWbox_metric(temperature, mesh_mask, depth_box=3000).item(),
)
print(
    "Deep water temperature:",
    temperature_DWbox_metric(temperature, mesh_mask, depth_box=3000).item(),
)

This returns:

Bottom water temperature: 4.666
Deep water temperature: 3.624