Four Box Model Use Tutorial

• ACTM Performer: JHU-APL;

• Author: Chace Ashcraft (chace.ashcraft@jhuapl.edu)

PyBAMOCS Tutorial

PyBAMOCS uses a simplified box model to simulate the Atlantic meridional overturning circulation (AMOC). The implementation is based on the paper by Gnanadesikan et al. (https://journals.ametsoc.org/view/journals/clim/31/22/jcli-d-18-0388.1.xml), which examined the stability of the AMOC. This tutorial assumes that you have downloaded the PyBAMOCs code from the PACMANs repository on GitHub and installed PyBAMOCS into your Python environment. Installation instructions for the PyBAMOCS module is provided below. For additional information about installing PyBAMOCS and the PyBAMOCS API, please see the PACMANS repository (https://github.com/JHUAPL/PACMANs), and navigate to box_model/python/.

This tutorial is meant to be a tool for getting started using the PyBAMOCS package and the included box_model functionality, and is not meant to be a comprehensive overview of all functions or capabilities of the included code. For a more thorough understanding, please consider viewing the code itself, and the examples in the scripts folder. There is also a more extensive tutorial in the form of a Jupyter Notebook in the notebooks directory of the PyBAMOCS module.

# PyBAMOCS Installation

We recommend installing PyBAMOCS in its own environment using venv (https://docs.python.org/3/library/venv.html) or conda (https://conda.io/projects/conda/en/latest/user-guide/install/index.html). Once in your new environment (if desired), perform the following steps to install the module.

• Clone the PACMANS repository (https://github.com/JHUAPL/PACMANs)

  • git clone https://github.com/JHUAPL/PACMANs.git

• Change directories to box_model/python/

  • cd box_model/python/

• Install locally using pip

  • pip install -e .

To execute the snippet:

Step 1: Follow the instructions above to install the PyBAMOCS module.

Step 2: The snippet itself may be run from the command line as follows:

` % python pybamocs_tutorial.py `

We encourage the interested user to examine the code below to identify the relevant objects and methods required to implement their specific use case.

# .............................................
# IMPORT STATEMENTS
# .............................................

import time
from matplotlib import pyplot as plt
from typing import Union

# Import the `box_model` function and its arguments
from pybamocs.box_model import box_model
from pybamocs.box_model import NORTH_IDX, SOUTH_IDX, LOW_IDX, DEEP_IDX
from pybamocs.box_model_args import (
    BoxModelBoxDimensions,
    BoxModelInitConditions,
    BoxModelTimeStep,
    BoxModelParameters,
)


# Convenience function
def print_box_model_argument_settings(argument: Union[BoxModelBoxDimensions, BoxModelInitConditions,
                                                      BoxModelTimeStep, BoxModelParameters]) -> None:
    """
    Print box model settings to console for a given box_model argument.
    :param argument: One of the box_model argument objects.
    """
    settings_dict = argument.to_dict()  # Get the settings as a python dictionary

    print("vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv")
    # Print each key, value pair in the dictionary
    for key, value in settings_dict.items():
        print(f"{key}={value}")
    print("^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^")
    print()


def main():
    print("Box model settings are divided into five groups:")
    # Each object represents a set of parameters to the model
    box_model_dimensions = BoxModelBoxDimensions()
    box_model_initial_conditions = BoxModelInitConditions()
    box_model_time_settings = BoxModelTimeStep()
    box_model_parameters = BoxModelParameters()

    # And there is one object for each group as shown above.
    print("------------------------------------------------------------")
    print("Box Dimensions:")
    print_box_model_argument_settings(box_model_dimensions)

    print("------------------------------------------------------------")
    print("Initial Conditions")
    print_box_model_argument_settings(box_model_initial_conditions)

    print("------------------------------------------------------------")
    print("Model Time Step Parameters:")
    print_box_model_argument_settings(box_model_time_settings)

    print("------------------------------------------------------------")
    print("Model Parameters")
    print_box_model_argument_settings(box_model_parameters)

    # To run the box model, simply call the `box_model` function with your argument objects
    print("Running the box model...")
    start = time.time()
    results = box_model(box_model_dimensions,
                        box_model_initial_conditions,
                        box_model_parameters,
                        box_model_time_settings)
    print(f"Run complete. Time: {time.time() - start} seconds.")

    # Note that the results are contained in a BoxModelResult object
    print("Box model output object:")
    print(results)
    print()

    # You can access specific results directly from the object
    print("Result details:")
    print("------------------------------------------------------")
    print("Shape of M_n for run:", results.M_n.shape)
    print()

    # or you can "unpack" them all at once
    M_n, M_upw, M_eddy, D_low, T, S, sigma0 = results.unpack()

    print("------------------------------------------------------")
    print("First 3 values of each output variable:")
    print("M_n:", M_n[:3])
    print("M_upw:", M_upw[:3])
    print("M_eddy:", M_eddy[:3])
    print("D_low:", D_low[:3])
    print("T_north:", T[NORTH_IDX, :3])
    print("T_south:", T[SOUTH_IDX, :3])
    print("T_low:", T[LOW_IDX, :3])
    print("T_deep:", T[DEEP_IDX, :3])
    print("S_north:", S[NORTH_IDX, :3])
    print("S_south:", S[SOUTH_IDX, :3])
    print("S_low:", S[LOW_IDX, :3])
    print("S_deep:", S[DEEP_IDX, :3])
    print("sigma0_north:", sigma0[NORTH_IDX, :3])
    print("sigma0_south:", sigma0[SOUTH_IDX, :3])
    print("sigma0_low:", sigma0[LOW_IDX, :3])
    print("sigma0_deep:", sigma0[DEEP_IDX, :3])

    print("Plotting results...")
    plt.plot(M_n, label='M_n')
    plt.plot(M_upw, label='M_upw')
    plt.plot(M_eddy, label='M_eddy')
    plt.plot(D_low, label='Dlow')
    plt.legend()
    plt.title("1D Box Model Outputs")

    fig, ax = plt.subplots(nrows=2, ncols=2)
    ax[0, 0].plot(T[NORTH_IDX], label='North')
    ax[0, 1].plot(T[SOUTH_IDX], label='South')
    ax[1, 0].plot(T[LOW_IDX], label='Low')
    ax[1, 1].plot(T[DEEP_IDX], label='Deep')
    ax[0, 0].legend()
    ax[0, 1].legend()
    ax[1, 0].legend()
    ax[1, 1].legend()
    plt.suptitle("Temperature in the 4 different boxes")
    plt.tight_layout()

    fig1, ax1 = plt.subplots(nrows=2, ncols=2)
    ax1[0, 0].plot(S[NORTH_IDX], label='North')
    ax1[0, 1].plot(S[SOUTH_IDX], label='South')
    ax1[1, 0].plot(S[LOW_IDX], label='Low')
    ax1[1, 1].plot(S[DEEP_IDX], label='Deep')
    ax1[0, 0].legend()
    ax1[0, 1].legend()
    ax1[1, 0].legend()
    ax1[1, 1].legend()
    plt.suptitle("Salinity in the 4 different boxes")
    plt.tight_layout()

    fig2, ax2 = plt.subplots(nrows=2, ncols=2)
    ax2[0, 0].plot(sigma0[NORTH_IDX], label='North')
    ax2[0, 1].plot(sigma0[SOUTH_IDX], label='South')
    ax2[1, 0].plot(sigma0[LOW_IDX], label='Low')
    ax2[1, 1].plot(sigma0[DEEP_IDX], label='Deep')
    ax2[0, 0].legend()
    ax2[0, 1].legend()
    ax2[1, 0].legend()
    ax2[1, 1].legend()
    plt.suptitle("Density in the 4 different boxes")
    plt.tight_layout()
    plt.show(block=False)

    # Consider an example of some data collection
    Fwn_values_to_test = [10000, 50000, 100000, 500000, 1000000]
    alternate_north_starting_temp = 4.0

    box_dims = BoxModelBoxDimensions()
    params = BoxModelParameters()
    time_step = BoxModelTimeStep()
    init = BoxModelInitConditions(T_north0=alternate_north_starting_temp)

    # Collect the data... (should take a few seconds)
    print(f"Generating data for alternate parameter settings: T_north={alternate_north_starting_temp}")
    start_time = time.time()
    results = []
    for fwn in Fwn_values_to_test:
        print(f"Fwn={fwn}")
        params.Fwn = fwn
        results.append(box_model(box_dims, init, params, time_step))
    time_to_collect_data = time.time() - start_time
    print(f"Total time to collect data: {time_to_collect_data} seconds")

    # Let's look at the difference in M_n for each run
    print("Plotting results...")
    plt.figure()
    for i in range(len(Fwn_values_to_test)):
        plt.plot(results[i].M_n, label=f"Fwn={Fwn_values_to_test[i]}")
    plt.title(f"M_n for different Northern Fluxes and T_north0={alternate_north_starting_temp}")
    plt.legend()
    plt.show()  # suspends script until plot windows are closed
    print("Done")


if __name__ == "__main__":
    main()