New commit following 1st revisions for GRL

data_for_figures/sam.csv

@@ -0,0 +1,42 @@
+Year,ERA-5,JRA55
+1979,1.05089,1.06014
+1980,-3.60272,-4.07698
+1981,-1.46379,-1.21394
+1982,-0.244872,-0.335743
+1983,0.535852,0.505942
+1984,-2.18357,-1.74108
+1985,1.15817,1.35848
+1986,-1.25842,-1.35748
+1987,-0.187979,-0.480559
+1988,-1.15005,-1.0362
+1989,2.07894,2.16804
+1990,-0.537585,-0.31202
+1991,-2.67013,-2.05847
+1992,-2.68876,-2.42423
+1993,1.31632,1.56319
+1994,-1.36703,-0.980318
+1995,-0.388389,0.0365294
+1996,-1.03008,-1.2204
+1997,-0.527492,-0.52408
+1998,4.6975,4.69334
+1999,2.71789,2.59908
+2000,-1.24927,-1.07009
+2001,1.44504,1.66743
+2002,-3.39209,-3.03674
+2003,-1.03079,-0.796827
+2004,0.397059,0.605366
+2005,-1.82684,-1.12621
+2006,0.705752,0.824439
+2007,-1.58114,-1.41207
+2008,1.72814,1.82205
+2009,-1.79592,-1.75158
+2010,3.47325,3.27941
+2011,-0.531007,-0.673399
+2012,0.850652,0.586529
+2013,0.257614,-0.119235
+2014,0.291953,-0.264715
+2015,3.1213,2.43896
+2016,1.66001,1.20798
+2017,1.92884,1.38202
+2018,1.5826,1.0196
+2019,-0.289823,-0.806082

final_scripts/sol_cfcs_sf6.jnl

@@ -0,0 +1,58 @@
+
+! define temperature and salinity variables
+let ztemp = ($1)
+let t_rk = (ztemp+273.16)*0.01
+let zsalt = ($2)
+
+
+!!! CFC-11 !!!
+
+! constants
+let A1 = -229.9261
+let A2 =  319.6552
+let A3 =  119.4471
+let A4 =  -1.39165
+let B1 =  -0.142382
+let B2 =   0.091459
+let B3 =  -0.0157274
+
+! solubility calculation
+let sol_cfc11_a = B1 + t_rk * (B2 + t_rk*B3)
+let sol_cfc11_b = exp(A1 + A2/t_rk + A3*ln(t_rk) + A4*t_rk*t_rk + zsalt * sol_cfc11_a)
+let sol_cfc11 = sol_cfc11_b * 1e3 * 1e-12  ! mol/l/atm --> mol/m3/pptv
+
+
+!!! CFC-12 !!!
+
+! constants
+let C1 = -218.0971 
+let C2 = 298.9702
+let C3 = 113.8049
+let C4 = -1.39165
+let D1 = -0.143566
+let D2 = 0.091015
+let D3 = -0.0153924
+
+! solubility calculation
+let sol_cfc12_a = D1 + t_rk * (D2 + t_rk*D3)
+let sol_cfc12_b = exp(C1 + C2/t_rk + C3*ln(t_rk) + C4*t_rk*t_rk + zsalt * sol_cfc12_a)
+let sol_cfc12 = sol_cfc12_b * 1e3 * 1e-12  ! mol/l/atm --> mol/m3/pptv
+
+
+!!! SF6 !!!
+
+! constants
+let E1 = -80.0343 
+let E2 = 117.232
+let E3 = 29.5817
+let E4 = 0.0
+let F1 = 0.0335183
+let F2 = -0.0373942
+let F3 = 0.00774862
+
+! solubility calculation
+let sol_sf6_a = F1 + t_rk * (F2 + t_rk*F3)
+let sol_sf6_b = exp(E1 + E2/t_rk + E3*ln(t_rk) + E4*t_rk*t_rk + zsalt * sol_sf6_a)
+let sol_sf6 = sol_sf6_b * 1e3 * 1e-12  ! mol/l/atm --> mol/m3/pptv
+
+

final_scripts/suppfig1.py

@@ -1,190 +1,137 @@
 # -*- coding: utf-8 -*-
 """
-Created on Thu Jan 21 13:45:34 2021
+Created on Tue Mar  2 10:39:44 2021
 
 Purpose
 -------
-    Figure of linear oxygen, solubility and AOU trends averaged between 200-1000 metres
-    in the hindcast model simulation and in the Helm climatology from 1970-1992
-    
-    
+    Create figure of change in the b-value of the Martin curve, which has been 
+    derived from the vertical profile of POC
+
+
 @author: pearseb
 """
 
 #%% imports
 
-from __future__ import unicode_literals
-
 import os
 import numpy as np
 import netCDF4 as nc
 import matplotlib.pyplot as plt
-import cmocean
-import cmocean.cm as cmo
+from matplotlib.gridspec import GridSpec
 import seaborn as sb
 sb.set(style='ticks')
 
+import cmocean
+import cmocean.cm as cmo
+
 from matplotlib.gridspec import GridSpec
 import matplotlib.ticker as mticker
 from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
 import cartopy
 import cartopy.crs as ccrs
 
 
-
 #%% get data
 
-os.chdir("C://Users//pearseb/Dropbox//PostDoc//my articles//historical model-data deoxygenation//data_for_figures")
-
-data = nc.Dataset("figure_compare_ito_jra55_ndep.nc", 'r')
-oxy_dep_dec = np.squeeze(data.variables['OXY_M_DEC'][...])
-tou_dep_dec = np.squeeze(data.variables['TOU_M_DEC'][...])*(-1)
+os.chdir("C://Users/pearseb/Dropbox/PostDoc/my articles/historical model-data deoxygenation/data_for_figures")
 
-data = nc.Dataset("figure_compare_ito_jra55.nc", 'r')
-oxy_pic_dec = np.squeeze(data.variables['OXY_M_DEC'][...])
-tou_pic_dec = np.squeeze(data.variables['TOU_M_DEC'][...])*(-1)
+data = nc.Dataset("map_hindcast_cmip6_ito_trends.nc",'r')
+oxy_ito = data.variables['OXY_ITO'][...]
+oxy_jra = data.variables['OXY_JRA55'][...]
+oxy_cmip = data.variables['OXY_CMIP'][...]
+oxy_cmip_agree = data.variables['OXY_CMIP_AGREE'][...]
 
 lon = data.variables['ETOPO60X'][...]
 lat = data.variables['ETOPO60Y'][...]
 
-data = nc.Dataset('figure_ndep_important.nc', 'r')
-ndep_important = data.variables['NDEP_IMPORTANT'][...]
-
-data = nc.Dataset("figure_compare_ito_jra55_6thpanel.nc", 'r')
-sol_trop = data.variables['SOL_TROP'][3:43]
-sol_extra = data.variables['SOL_EXTRA'][3:43]
-sol_glob = data.variables['SOL_GLOB'][3:43]
-oxy_pic_trop = data.variables['OXY_PIC_TROP'][3:43]
-oxy_pic_extra = data.variables['OXY_PIC_EXTRA'][3:43]
-oxy_pic_glob = data.variables['OXY_PIC_GLOB'][3:43]
-aou_pic_trop = data.variables['AOU_PIC_TROP'][3:43]*(-1)
-aou_pic_extra = data.variables['AOU_PIC_EXTRA'][3:43]*(-1)
-aou_pic_glob = data.variables['AOU_PIC_GLOB'][3:43]*(-1)
-oxy_dep_trop = data.variables['OXY_DEP_TROP'][3:43]
-oxy_dep_extra = data.variables['OXY_DEP_EXTRA'][3:43]
-oxy_dep_glob = data.variables['OXY_DEP_GLOB'][3:43]
-aou_dep_trop = data.variables['AOU_DEP_TROP'][3:43]*(-1)
-aou_dep_extra = data.variables['AOU_DEP_EXTRA'][3:43]*(-1)
-aou_dep_glob = data.variables['AOU_DEP_GLOB'][3:43]*(-1)
-
-years = np.arange(1975,2015,1)
 
 data.close()
 
 
-#%% figure specifics 
+
+#%% figure specifics
 
 fslab = 15
 fstic = 13
 lw = 0.75
 gridalf = 0.5
 
 colmap1 = cmocean.tools.lighten(cmo.balance, 0.8)
-levs1 = np.array([-0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5])
-cont1 = levs1[::2]
+levs1 = np.arange(-4, 4.1, 0.25)
+
 
 contcol = 'black'
-contwid = 0.5
+contwid=0.75
 
 proj = ccrs.Robinson(central_longitude=200)
 lons,lats = np.meshgrid(lon,lat)
 
 
-sol_cont = np.arange(-50,51,1)
-aou_cont = np.arange(-50,51,1)
-sol_cont, aou_cont = np.meshgrid(sol_cont, aou_cont)
-oxy_cont = sol_cont+aou_cont
-
 
-#%% make figure
+#%% figure
 
-fig = plt.figure(figsize=(6,9))
-gs = GridSpec(15,1)
+fig = plt.figure(figsize=(7,9.5))
+gs = GridSpec(3,1)
 
-ax1 = plt.subplot(gs[0:6], projection=proj)
+ax1 = plt.subplot(gs[0,0], projection=proj)
 ax1.tick_params(labelsize=fstic)
 ax1.gridlines(linestyle='--', linewidth=0.5, color='grey', alpha=gridalf, zorder=3)
 ax1.add_feature(cartopy.feature.LAND, zorder=1, facecolor='silver')
 ax1.coastlines()
-p1 = plt.contourf(lons, lats, oxy_dep_dec-oxy_pic_dec, transform=ccrs.PlateCarree(), cmap=colmap1, levels=levs1, vmin=np.min(levs1),vmax=np.max(levs1), zorder=1, extend='both')
-plt.contourf(lons,lats, ndep_important, transform=ccrs.PlateCarree(), colors='none', levels=[0.9,1.1], hatches=['......'],zorder=2)
-c1 = plt.contour(lons, lats, oxy_dep_dec-oxy_pic_dec, transform=ccrs.PlateCarree(), colors=contcol, linewidths=contwid, levels=cont1, zorder=2)
-
-aou_dep_change = aou_dep_glob - aou_dep_glob[0]
-aou_pic_change = aou_pic_glob - aou_pic_glob[0]
-sol_change =  sol_glob - sol_glob[0]
-
-ax2 = plt.subplot(gs[6:9,0])
-ax2.tick_params(labelsize=fstic, labelbottom=False)
-plt.plot(aou_pic_change, sol_change, color='k', label='Preindustrial N$_r$ deposition')
-plt.plot(aou_dep_change, sol_change, color='forestgreen', label='Anthropogenic N$_r$ deposition')
-c1 = plt.contour(aou_cont, sol_cont, oxy_cont, levels=np.arange(-100,101,0.2), colors=contcol, linewidths=contwid)
-plt.ylim(-0.2,0.2)
-plt.xlim(-1.5,0.05)
-plt.ylabel('$\Delta$ O$_2^{sol}$ ($\mu$M)', fontsize=fslab)
-plt.legend(frameon=False, loc='upper left', ncol=1, labelspacing=0.05, bbox_to_anchor=(0.2,1.45))
-
-aou_dep_change = aou_dep_trop - aou_dep_trop[0]
-aou_pic_change = aou_pic_trop - aou_pic_trop[0]
-sol_change =  sol_trop - sol_trop[0]
-
-ax3 = plt.subplot(gs[9:12,0])
-ax3.tick_params(labelsize=fstic, labelbottom=False)
-plt.plot(aou_pic_change, sol_change, color='k', label='Preindustrial N$_r$ deposition')
-plt.plot(aou_dep_change, sol_change, color='forestgreen', label='Anthropogenic N$_r$ deposition')
-c2 = plt.contour(aou_cont, sol_cont, oxy_cont, levels=np.arange(-100,101,0.2), colors=contcol, linewidths=contwid)
-plt.ylim(-0.25,0.25)
-plt.xlim(-1.5,0.05)
-plt.ylabel('$\Delta$ O$_2^{sol}$ ($\mu$M)', fontsize=fslab)
-
-aou_dep_change = aou_dep_extra - aou_dep_extra[0]
-aou_pic_change = aou_pic_extra - aou_pic_extra[0]
-sol_change =  sol_extra - sol_extra[0]
-
-ax4 = plt.subplot(gs[12:15,0])
-ax4.tick_params(labelsize=fstic, labelbottom=True)
-plt.plot(aou_pic_change, sol_change, color='k', label='Preindustrial N$_r$ deposition')
-plt.plot(aou_dep_change, sol_change, color='forestgreen', label='Anthropogenic N$_r$ deposition')
-c3 = plt.contour(aou_cont, sol_cont, oxy_cont, levels=np.arange(-100,101,0.2), colors=contcol, linewidths=contwid)
-plt.ylim(-0.25,0.25)
-plt.xlim(-1.5,0.05)
-plt.ylabel('$\Delta$ O$_2^{sol}$ ($\mu$M)', fontsize=fslab)
-plt.xlabel('$\Delta$ O$_2^{AOU}$ ($\mu$M)', fontsize=fslab)
-
-
-fig.subplots_adjust(top=0.875, bottom=0.1, left=0.15, right=0.9, hspace=1.5)
-
-
-xx = 0.5; yy = 1.05
-plt.text(xx, yy, '$\Delta$ O$_2$ due to N$_r$ deposition', fontsize=fstic, transform=ax1.transAxes, va='center', ha='center', fontweight='bold')
-
-xx = 1.05; yy = 0.5
-plt.text(xx,yy,'Global', transform=ax2.transAxes, rotation=90, ha='center', va='center', fontsize=fstic)
-plt.text(xx,yy,'20$^{\circ}$S-20$^{\circ}$N', transform=ax3.transAxes, rotation=90, ha='center', va='center', fontsize=fstic)
-plt.text(xx,yy,'Extra-tropics', transform=ax4.transAxes, rotation=90, ha='center', va='center', fontsize=fstic)
+p1 = plt.contourf(lons, lats, oxy_ito, transform=ccrs.PlateCarree(), cmap=colmap1, levels=levs1, vmin=np.min(levs1),vmax=np.max(levs1), zorder=1, extend='both')
+c1 = plt.contour(lons, lats, oxy_ito, transform=ccrs.PlateCarree(), levels=[0], zorder=2, colors=contcol, linewidths=contwid)
+
+ax2 = plt.subplot(gs[1,0], projection=proj)
+ax2.tick_params(labelsize=fstic)
+ax2.gridlines(linestyle='--', linewidth=0.5, color='grey', alpha=gridalf, zorder=3)
+ax2.add_feature(cartopy.feature.LAND, zorder=1, facecolor='silver')
+ax2.coastlines()
+p2 = plt.contourf(lons, lats, oxy_jra, transform=ccrs.PlateCarree(), cmap=colmap1, levels=levs1, vmin=np.min(levs1),vmax=np.max(levs1), zorder=1, extend='both')
+c2 = plt.contour(lons, lats, oxy_jra, transform=ccrs.PlateCarree(), levels=[0], zorder=2, colors=contcol, linewidths=contwid)
+
+ax3 = plt.subplot(gs[2,0], projection=proj)
+ax3.tick_params(labelsize=fstic)
+ax3.gridlines(linestyle='--', linewidth=0.5, color='grey', alpha=gridalf, zorder=3)
+ax3.add_feature(cartopy.feature.LAND, zorder=1, facecolor='silver')
+ax3.coastlines()
+p3 = plt.contourf(lons, lats, oxy_cmip, transform=ccrs.PlateCarree(), cmap=colmap1, levels=levs1, vmin=np.min(levs1),vmax=np.max(levs1), zorder=2, extend='both')
+c3 = plt.contour(lons, lats, oxy_cmip, transform=ccrs.PlateCarree(), levels=[0], zorder=2, colors=contcol, linewidths=contwid)
+plt.contourf(lons, lats, oxy_cmip_agree, transform=ccrs.PlateCarree(), colors='none', levels=[0,0.5,1.5], hatches=[' ','....'], zorder=1)
+
+fig.subplots_adjust(top=0.95, bottom=0.05, left=0.05, right=0.82, wspace=0.05, hspace=0.05)
+
+
+#xx = 0.5; yy = 1.05
+#plt.text(xx, yy, 'Ideal age', fontsize=fslab, transform=ax1.transAxes, va='center', ha='center', fontweight='bold')
+#plt.text(xx, yy, 'Biological demand', fontsize=fslab, transform=ax2.transAxes, va='center', ha='center', fontweight='bold')
 
 xx = 0.05; yy = 0.95
 plt.text(xx, yy, 'a', fontsize=fslab+2, transform=ax1.transAxes, va='center', ha='center', fontweight='bold')
-plt.text(xx-0.025, yy+0.15, 'b', fontsize=fslab+2, transform=ax2.transAxes, va='center', ha='center', fontweight='bold')
-plt.text(xx-0.025, yy+0.15, 'c', fontsize=fslab+2, transform=ax3.transAxes, va='center', ha='center', fontweight='bold')
-plt.text(xx-0.025, yy+0.15, 'd', fontsize=fslab+2, transform=ax4.transAxes, va='center', ha='center', fontweight='bold')
+plt.text(xx, yy, 'b', fontsize=fslab+2, transform=ax2.transAxes, va='center', ha='center', fontweight='bold')
+plt.text(xx, yy, 'c', fontsize=fslab+2, transform=ax3.transAxes, va='center', ha='center', fontweight='bold')
 
 
-cbax1 = fig.add_axes([0.2, 0.9, 0.65, 0.025])
-cbar1 = plt.colorbar(p1, cax=cbax1, orientation='horizontal', ticks=levs1[::2])
-cbax1.set_xlabel('$\Delta$ O$_2$ ($\mu$M decade$^{-1}$)', fontsize=fslab)
-cbax1.tick_params(labelsize=fstic, top=True, labeltop=True, bottom=False, labelbottom=False)
-cbax1.xaxis.set_label_position('top')
+cbax1 = fig.add_axes([0.84, 0.68, 0.03, 0.25])
+cbar1 = plt.colorbar(p1, cax=cbax1, orientation='vertical', ticks=levs1[::4])
+cbax1.tick_params(labelsize=fstic)
+cbax1.set_ylabel("$\Delta$ O$_2$ ($\mu$M decade$^{-1}$)", fontsize=fslab)
 
+cbax2 = fig.add_axes([0.84, 0.375, 0.03, 0.25])
+cbar2 = plt.colorbar(p2, cax=cbax2, orientation='vertical', ticks=levs1[::4])
+cbax2.tick_params(labelsize=fstic)
+cbax2.set_ylabel("$\Delta$ O$_2$ ($\mu$M decade$^{-1}$)", fontsize=fslab)
 
-#%% save figure
+cbax3 = fig.add_axes([0.84, 0.07, 0.03, 0.25])
+cbar3 = plt.colorbar(p3, cax=cbax3, orientation='vertical', ticks=levs1[::4])
+cbax3.tick_params(labelsize=fstic)
+cbax3.set_ylabel("$\Delta$ O$_2$ ($\mu$M decade$^{-1}$)", fontsize=fslab)
 
-plt.clabel(c1, manual=True, fmt='%.1f', fontsize=fstic-2)
-plt.clabel(c2, manual=True, fmt='%.1f', fontsize=fstic-2)
-plt.clabel(c3, manual=True, fmt='%.1f', fontsize=fstic-2)
 
 
+#%% save figure
+
 os.chdir("C://Users//pearseb/Dropbox//PostDoc//my articles//historical model-data deoxygenation//final_figures")
 fig.savefig('fig-suppfig1.png', dpi=300, bbox_inches='tight')
 fig.savefig('fig-suppfig1_trans.png', dpi=300, bbox_inches='tight', transparent=True)
 
+