{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Code for Fig. 3 and fig. S3 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", "dir_hr = '/glade/work/yeager/iHesp/B1850/'\n", "dir_lr = '/glade/work/yeager/iHesp/B1850_LR/'\n", "f1_hr = f'{datadir}/B.E.13.B1850C5.ne120_t12.sehires38.003.sunway_02.pop.h.015001_050012.MOCsig.osnap.nc'\n", "f1_lr = f'{datadir}/B.E.13.B1850C5.ne30g16.sehires38.003.sunway.pop.h.000101_050012.MOCsig.osnap.nc'\n", "f2_hr = f'{datadir}/B.E.13.B1850C5.ne120_t12.sehires38.003.sunway_02.pop.h.020001_050012.MOCsig.nc'\n", "f2_lr = f'{datadir}/B.E.13.B1850C5.ne30g16.sehires38.003.sunway.pop.h.000101_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", "f3_lr = f'{datadir}/B.E.13.B1850C5.ne30g16.sehires38.003.sunway.pop.h.000101_050012.WMF.nc'\n", "\n", "ds1_hr = open_POPdataset(f1_hr) \n", "ds1_lr = open_POPdataset(f1_lr)\n", "ds2_hr = open_POPdataset(f2_hr) \n", "ds2_lr = open_POPdataset(f2_lr)\n", "ds3_hr = open_POPdataset(f3_hr) \n", "ds3_lr = open_POPdataset(f3_lr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute Anomalies and Std Dev (years 200-500)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "ds1_hr_ann = ds1_hr.groupby('time.year').mean('time')\n", "ds1_lr_ann = ds1_lr.groupby('time.year').mean('time')\n", "ds2_hr_ann = ds2_hr.groupby('time.year').mean('time')\n", "ds2_lr_ann = ds2_lr.groupby('time.year').mean('time')\n", "ds3_hr_ann = ds3_hr.groupby('time.year').mean('time')\n", "ds3_lr_ann = ds3_lr.groupby('time.year').mean('time')\n", "\n", "ds1_hr_ann=ds1_hr_ann.rename({'year':'time'}).sel(time=slice(200,501))\n", "ds1_lr_ann=ds1_lr_ann.rename({'year':'time'}).sel(time=slice(200,501))\n", "ds2_hr_ann=ds2_hr_ann.rename({'year':'time'}).sel(time=slice(200,501))\n", "ds2_lr_ann=ds2_lr_ann.rename({'year':'time'}).sel(time=slice(200,501))\n", "ds3_hr_ann=ds3_hr_ann.rename({'year':'time'}).sel(time=slice(200,501))\n", "ds3_lr_ann=ds3_lr_ann.rename({'year':'time'}).sel(time=slice(200,501))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Annual Anomalies\n", "ds1_lr_mean = ds1_lr_ann.mean('time')\n", "ds1_lr_annanom = ds1_lr_ann - ds1_lr_mean\n", "ds1_hr_mean = ds1_hr_ann.mean('time')\n", "ds1_hr_annanom = ds1_hr_ann - ds1_hr_mean\n", "ds2_lr_mean = ds2_lr_ann.mean('time')\n", "ds2_lr_annanom = ds2_lr_ann - ds2_lr_mean\n", "ds2_hr_mean = ds2_hr_ann.mean('time')\n", "ds2_hr_annanom = ds2_hr_ann - ds2_hr_mean\n", "ds3_lr_mean = ds3_lr_ann.mean('time')\n", "ds3_lr_annanom = ds3_lr_ann - ds3_lr_mean\n", "ds3_hr_mean = ds3_hr_ann.mean('time')\n", "ds3_hr_annanom = ds3_hr_ann - ds3_hr_mean\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Annual Detrended Anomalies\n", "ds1_lr_annanom_dt = xr.apply_ufunc(signal.detrend, ds1_lr_annanom.fillna(0), kwargs={'axis': 0}).where(ds1_lr_annanom.notnull())\n", "ds2_lr_annanom_dt = xr.apply_ufunc(signal.detrend, ds2_lr_annanom.fillna(0), kwargs={'axis': 0}).where(ds2_lr_annanom.notnull())\n", "ds3_lr_annanom_dt = xr.apply_ufunc(signal.detrend, ds3_lr_annanom.fillna(0), kwargs={'axis': 0}).where(ds3_lr_annanom.notnull())\n", "\n", "ds1_hr_annanom_dt = xr.apply_ufunc(signal.detrend, ds1_hr_annanom.fillna(0), kwargs={'axis': 0}).where(ds1_hr_annanom.notnull())\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())\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Temporal Filtering" ] }, { "cell_type": "code", "execution_count": 7, "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": 8, "metadata": {}, "outputs": [], "source": [ "# Low-Pass (lp) anomalies:\n", "ds1_lr_lpanom_dt = xr.apply_ufunc(signal.sosfiltfilt, filtsos, ds1_lr_annanom_dt.fillna(0), kwargs={'padtype':'even','axis':0}).where(ds1_lr_annanom_dt.notnull())\n", "ds2_lr_lpanom_dt = xr.apply_ufunc(signal.sosfiltfilt, filtsos, ds2_lr_annanom_dt.fillna(0), kwargs={'padtype':'even','axis':0}).where(ds2_lr_annanom_dt.notnull())\n", "ds3_lr_lpanom_dt = xr.apply_ufunc(signal.sosfiltfilt, filtsos, ds3_lr_annanom_dt.fillna(0), kwargs={'padtype':'even','axis':0}).where(ds3_lr_annanom_dt.notnull())\n", "\n", "ds1_hr_lpanom_dt = xr.apply_ufunc(signal.sosfiltfilt, filtsos, ds1_hr_annanom_dt.fillna(0), kwargs={'padtype':'even','axis':0}).where(ds1_hr_annanom_dt.notnull())\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())\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Detrended, filtered Std Dev:\n", "ds1_lr_anndtstd = ds1_lr_annanom_dt.std('time')\n", "ds1_hr_anndtstd = ds1_hr_annanom_dt.std('time')\n", "ds1_lr_lpdtstd = ds1_lr_lpanom_dt.std('time')\n", "ds1_hr_lpdtstd = ds1_hr_lpanom_dt.std('time')\n", "\n", "ds2_lr_anndtstd = ds2_lr_annanom_dt.std('time')\n", "ds2_hr_anndtstd = ds2_hr_annanom_dt.std('time')\n", "ds2_lr_lpdtstd = ds2_lr_lpanom_dt.std('time')\n", "ds2_hr_lpdtstd = ds2_hr_lpanom_dt.std('time')\n", "\n", "ds3_lr_anndtstd = ds3_lr_annanom_dt.std('time')\n", "ds3_hr_anndtstd = ds3_hr_annanom_dt.std('time')\n", "ds3_lr_lpdtstd = ds3_lr_lpanom_dt.std('time')\n", "ds3_hr_lpdtstd = ds3_hr_lpanom_dt.std('time')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Fig 3" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
<xarray.DataArray 'wmf_region' (wmf_region: 7)>\n", "array(['All (>0)', 'SPG_west (1)', 'SPG_east (2)', 'Irminger Sea (3)',\n", " 'Labrador Sea (4)', 'Norwegian Sea (5)', 'Arctic (6)'], dtype=object)\n", "Coordinates:\n", " * wmf_region (wmf_region) object 'All (>0)' 'SPG_west (1)' ... 'Arctic (6)'
array(['All (>0)', 'SPG_west (1)', 'SPG_east (2)', 'Irminger Sea (3)',\n", " 'Labrador Sea (4)', 'Norwegian Sea (5)', 'Arctic (6)'], dtype=object)
array(['All (>0)', 'SPG_west (1)', 'SPG_east (2)', 'Irminger Sea (3)',\n", " 'Labrador Sea (4)', 'Norwegian Sea (5)', 'Arctic (6)'], dtype=object)