{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Code for Fig. 5 and fig. S5 of\n", "#### Yeager et al., 2021: An Outsized Role for the Labrador Sea in the Multidecadal Variability of the Atlantic Overturning Circulation, *Science Advances*." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import xarray as xr \n", "import numpy as np \n", "import cftime\n", "import copy\n", "import scipy.stats\n", "from scipy import signal\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def open_POPdataset(x):\n", " ds = xr.open_dataset(x,decode_times=False)\n", " attrs=ds.time.attrs.copy()\n", " ds = ds.assign_coords(time=ds.time.values - 15)\n", " ds.time.attrs = attrs\n", " ds = xr.decode_cf(ds)\n", " return ds" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "datadir = '/glade/scratch/yeager/YeagerEA_ScienceAdvances_2021/'\n", "f2_hr = f'{datadir}/B.E.13.B1850C5.ne120_t12.sehires38.003.sunway_02.pop.h.020001_050012.MOCsig.nc'\n", "f3_hr = f'{datadir}/B.E.13.B1850C5.ne120_t12.sehires38.003.sunway_02.pop.h.020001_050012.WMF.nc'\n", "\n", "ds2_hr = open_POPdataset(f2_hr) \n", "ds3_hr = open_POPdataset(f3_hr) " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Add additional regions to WMF datasets\n", "ds3_hr1 = ds3_hr.isel(wmf_region=[1,2,3,5,6]).sum('wmf_region')\n", "ds3_hr1 = ds3_hr1.assign_coords({'wmf_region':'ALL minus LAB (7)'})\n", "ds3_hr2 = ds3_hr.isel(wmf_region=[1,2,3]).sum('wmf_region')\n", "ds3_hr2 = ds3_hr2.assign_coords({'wmf_region':'IRM+SPG (8)'})\n", "ds3_hr3 = ds3_hr.isel(wmf_region=[1,4]).sum('wmf_region')\n", "ds3_hr3 = ds3_hr3.assign_coords({'wmf_region':'LAB+SPG-west (9)'})\n", "ds3_hr4 = ds3_hr.isel(wmf_region=[2,3]).sum('wmf_region')\n", "ds3_hr4 = ds3_hr4.assign_coords({'wmf_region':'IRM+SPG-east (10)'})\n", "ds3_hr = xr.concat([ds3_hr,ds3_hr1,ds3_hr2,ds3_hr3,ds3_hr4],dim='wmf_region')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute Anomalies and Std Dev (years 200-500)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "ds2_hr_ann = ds2_hr.groupby('time.year').mean('time')\n", "ds3_hr_ann = ds3_hr.groupby('time.year').mean('time')\n", "\n", "ds2_hr_ann=ds2_hr_ann.rename({'year':'time'}).sel(time=slice(200,501))\n", "ds3_hr_ann=ds3_hr_ann.rename({'year':'time'}).sel(time=slice(200,501))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Annual Anomalies\n", "ds2_hr_mean = ds2_hr_ann.mean('time')\n", "ds2_hr_annanom = ds2_hr_ann - ds2_hr_mean\n", "ds3_hr_mean = ds3_hr_ann.mean('time')\n", "ds3_hr_annanom = ds3_hr_ann - ds3_hr_mean" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Annual Detrended Anomalies\n", "ds2_hr_annanom_dt = xr.apply_ufunc(signal.detrend, ds2_hr_annanom.fillna(0), kwargs={'axis': 0}).where(ds2_hr_annanom.notnull())\n", "ds3_hr_annanom_dt = xr.apply_ufunc(signal.detrend, ds3_hr_annanom.fillna(0), kwargs={'axis': 0}).where(ds3_hr_annanom.notnull())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Temporal Filtering" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "cutoff= 10.0 years\n" ] } ], "source": [ "# 10-year butterworth low-pass filter\n", "fs=1/(365*24*3600) # 1 year in Hz (sampling frequency)\n", "nyquist = fs / 2 # 0.5 times the sampling frequency\n", "cutoff = fs/10 # 10-year cutoff frequency\n", "cutoff = cutoff/nyquist # as fraction of nyquist \n", "print('cutoff= ',(1/(cutoff*nyquist))/(365*24*3600),' years') \n", "filtsos = signal.butter(4, cutoff, 'lowpass', output='sos') #low pass filter\n", "filtb, filta = signal.butter(4, cutoff, 'lowpass')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Now Apply Filter Using Xarray apply_ufunc" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Low-Pass (lp) anomalies:\n", "ds2_hr_lpanom_dt = xr.apply_ufunc(signal.sosfiltfilt, filtsos, ds2_hr_annanom_dt.fillna(0), kwargs={'padtype':'even','axis':0}).where(ds2_hr_annanom_dt.notnull())\n", "ds3_hr_lpanom_dt = xr.apply_ufunc(signal.sosfiltfilt, filtsos, ds3_hr_annanom_dt.fillna(0), kwargs={'padtype':'even','axis':0}).where(ds3_hr_annanom_dt.notnull())" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Detrended, filtered Std Dev:\n", "ds2_hr_anndtstd = ds2_hr_annanom_dt.std('time')\n", "ds2_hr_lpdtstd = ds2_hr_lpanom_dt.std('time')\n", "\n", "ds3_hr_anndtstd = ds3_hr_annanom_dt.std('time')\n", "ds3_hr_lpdtstd = ds3_hr_lpanom_dt.std('time')" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# First, create low-pass-filtered, detrended AMOC that includes the mean\n", "ds2_hr_lpann_dt = ds2_hr_lpanom_dt + ds2_hr_mean\n", "ds3_hr_lpann_dt = ds3_hr_lpanom_dt + ds3_hr_mean" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### LSW density range determined here:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([-0.15392032, 0.37033764, 1.42853378, 1.93606328, 0.55137143,\n", " 0.01353341]),\n", " array([36.925, 36.975, 37.025, 37.075, 37.125, 37.175]))" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds3_hr_lpann_dt.isel(wmf_region=4).mean('time').WMF.sel(sigma_wmf=slice(36.9,37.2)).values,ds3_hr_lpann_dt.sel(sigma_wmf=slice(36.9,37.2)).sigma_wmf.values" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "hr_lsw = [36.95, 37.175]\n", "hr_dlsw = [37.0625, 37.175]" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
<xarray.DataArray 'wmf_region' (wmf_region: 11)>\n", "array(['All (>0)', 'SPG_west (1)', 'SPG_east (2)', 'Irminger Sea (3)',\n", " 'Labrador Sea (4)', 'Norwegian Sea (5)', 'Arctic (6)',\n", " 'ALL minus LAB (7)', 'IRM+SPG (8)', 'LAB+SPG-west (9)',\n", " 'IRM+SPG-east (10)'], dtype=object)\n", "Coordinates:\n", " * wmf_region (wmf_region) object 'All (>0)' ... 'IRM+SPG-east (10)'
array(['All (>0)', 'SPG_west (1)', 'SPG_east (2)', 'Irminger Sea (3)',\n", " 'Labrador Sea (4)', 'Norwegian Sea (5)', 'Arctic (6)',\n", " 'ALL minus LAB (7)', 'IRM+SPG (8)', 'LAB+SPG-west (9)',\n", " 'IRM+SPG-east (10)'], dtype=object)
array(['All (>0)', 'SPG_west (1)', 'SPG_east (2)', 'Irminger Sea (3)',\n", " 'Labrador Sea (4)', 'Norwegian Sea (5)', 'Arctic (6)',\n", " 'ALL minus LAB (7)', 'IRM+SPG (8)', 'LAB+SPG-west (9)',\n", " 'IRM+SPG-east (10)'], dtype=object)