From 0816e7dafc3b7d100e8db555d5ac8269d9792040 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Fri, 12 Feb 2016 13:09:30 -0500 Subject: [PATCH 1/9] Analysis script to compute y-direction heat transport Initially works with T_ady_2d, the advective component. The full heat transport should also inclide T_diffy_2d. If the script is invoked with only T_ady_2d, a warning message is printed to the screen and annoateed to the plot. --- tools/analysis/poleward_heat_transport.py | 117 ++++++++++++++++++++++ 1 file changed, 117 insertions(+) create mode 100755 tools/analysis/poleward_heat_transport.py diff --git a/tools/analysis/poleward_heat_transport.py b/tools/analysis/poleward_heat_transport.py new file mode 100755 index 0000000000..8a05209921 --- /dev/null +++ b/tools/analysis/poleward_heat_transport.py @@ -0,0 +1,117 @@ +#!/usr/bin/env python + +import netCDF4 +import numpy +import m6plot +import matplotlib.pyplot as plt +import sys +import warnings + +def run(): + try: import argparse + except: raise Exception('This version of python is not new enough. python 2.7 or newer is required.') + parser = argparse.ArgumentParser(description='''Script for plotting plotting poleward heat transport.''') + parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'T_ady_2d' and 'T_diffy_2d'.''') + parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') + parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') + parser.add_argument('-g','--gridspecdir', type=str, required=True, + help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') + cmdLineArgs = parser.parse_args() + main(cmdLineArgs) + +def main(cmdLineArgs,stream=None): + #x = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['x'][::2,::2] + y = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['y'][::2,::2] + msk = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_mask.nc').variables['mask'][:] + area = msk*netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_topog.nc').variables['depth'][:] + basin_code = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/basin_codes.v20140629.nc').variables['basin'][:] + + if stream != None: + if len(stream) != 2: + raise ValueError('If specifying output streams, exactly two streams are needed for this analysis') + + rootGroup = netCDF4.MFDataset( cmdLineArgs.annual_file ) + if 'T_ady_2d' in rootGroup.variables: + varName = 'T_ady_2d' + if len(rootGroup.variables[varName].shape)==4: advective = rootGroup.variables[varName][:].mean(axis=0).filled(0.) + else: advective = rootGroup.variables[varName][:].filled(0.) + else: raise Exception('Could not find "T_ady_2d" in file "%s"'%(cmdLineArgs.annual_file)) + if 'T_diffy_2d' in rootGroup.variables: + varName = 'T_diffy_2d' + if len(rootGroup.variables[varName].shape)==4: diffusive = rootGroup.variables[varName][:].mean(axis=0).filled(0.) + else: diffusive = rootGroup.variables[varName][:].filled(0.) + else: + diffusive = None + warnings.warn('Diffusive temperature term not found. This will result in an underestimation of the heat transport.') + + def heatTrans(advective, diffusive=None, vmask=None): + """Converts vertically integrated temperature advection into heat transport""" + rho0 = 1.035e3; Cp = 3989. + if diffusive != None: HT = advective + diffusive + else: HT = advective + HT = HT * (rho0 * Cp); HT = HT * 1.e-15 # convert to PW + if vmask != None: HT = HT*vmask + HT = HT.sum(axis=-1); HT = HT.squeeze() # sum in x-direction + return HT + + def plotHeatTrans(y, HT, title): + plt.plot(y, y*0., 'k', linewidth=0.5) + plt.plot(y, HT, 'k', linewidth=1.5) + plt.xlim(-80,90); plt.ylim(-2.0,3.0) + plt.title(title) + plt.grid(True) + + def annotatePlot(label): + fig = plt.gcf() + fig.text(0.1,0.85,label) + + m6plot.setFigureSize(npanels=1) + + # Global Heat Transport + HTplot = heatTrans(advective,diffusive) + yy = y[1:,:].max(axis=-1) + plotHeatTrans(yy,HTplot,title='Global Y-Direction Heat Transport [PW]') + plt.xlabel(r'Latitude [$\degree$N]') + plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) + if diffusive == None: annotatePlot('Warning: Diffusive component of transport is missing.') + if stream != None: + plt.savefig(stream[0]) + else: + plt.savefig(cmdLineArgs.outdir+'/HeatTransport_global.png') + plt.show(block=False) + + # Atlantic Heat Transport + plt.clf() + m = 0*basin_code; m[(basin_code==2) | (basin_code==4)] = 1 + HTplot = heatTrans(advective, diffusive, vmask=m*numpy.roll(m,1,axis=1)) + yy = y[1:,:].max(axis=-1) + HTplot[yy<-34] = numpy.nan + plotHeatTrans(yy,HTplot,title='Atlantic Y-Direction Heat Transport [PW]') + plt.xlabel(r'Latitude [$\degree$N]') + plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) + if diffusive == None: annotatePlot('Warning: Diffusive component of transport is missing.') + if stream != None: + plt.savefig(stream[0]) + else: + plt.savefig(cmdLineArgs.outdir+'/HeatTransport_Atlantic.png') + plt.show(block=False) + + # Indo-Pacific Heat Transport + plt.clf() + m = 0*basin_code; m[(basin_code==3) | (basin_code==5)] = 1 + HTplot = heatTrans(advective, diffusive, vmask=m*numpy.roll(m,1,axis=1)) + yy = y[1:,:].max(axis=-1) + HTplot[yy<-38] = numpy.nan + plotHeatTrans(yy,HTplot,title='Indo-Pacific Y-Direction Heat Transport [PW]') + plt.xlabel(r'Latitude [$\degree$N]') + if diffusive == None: annotatePlot('Warning: Diffusive component of transport is missing.') + plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) + if stream != None: + plt.savefig(stream[0]) + else: + plt.savefig(cmdLineArgs.outdir+'/HeatTransport_IndoPac.png') + plt.show(block=False) + +if __name__ == '__main__': + run() From c344903257329ff7f03277226521dfa74ecf3a0d Mon Sep 17 00:00:00 2001 From: John Krasting Date: Thu, 18 Feb 2016 08:17:01 -0500 Subject: [PATCH 2/9] Analysis mods for use with multiple grids Feature enhancements: - Ability to read either gridspec directory or tar file - On-line basin mask generation if basin mask file is missing Added a Makefile to run tests of the common frepp analysis scripts on both the 0.25 and 0.5 degree configurations. --- tools/analysis/MLD_003.py | 21 ++- tools/analysis/Makefile | 47 ++++++ tools/analysis/SSS_bias_WOA05.py | 29 +++- tools/analysis/SST_bias_WOA05.py | 30 +++- tools/analysis/basin_test.py | 54 ++++++ tools/analysis/depth_average_T_bias.py | 29 +++- tools/analysis/m6toolbox.py | 154 +++++++++++++++++- tools/analysis/meridional_overturning.py | 43 +++-- tools/analysis/poleward_heat_transport.py | 54 +++--- .../vertical_sections_annual_bias_WOA05.py | 33 +++- tools/analysis/zonal_S_bias_WOA05.py | 32 +++- tools/analysis/zonal_T_bias_WOA05.py | 29 +++- 12 files changed, 482 insertions(+), 73 deletions(-) create mode 100644 tools/analysis/Makefile create mode 100644 tools/analysis/basin_test.py diff --git a/tools/analysis/MLD_003.py b/tools/analysis/MLD_003.py index 272d6d85ab..31b1030e87 100755 --- a/tools/analysis/MLD_003.py +++ b/tools/analysis/MLD_003.py @@ -3,7 +3,9 @@ import netCDF4 import numpy import m6plot +import m6toolbox import matplotlib.pyplot as plt +import os try: import argparse except: raise Exception('This version of python is not new enough. python 2.7 or newer is required.') @@ -15,17 +17,26 @@ parser.add_argument('-od','--obsdata', type=str, default='/archive/gold/datasets/obs/Hosada2010_MLD_climatology.v20140515.nc', help='''File containing the observational MLD data (Hosoda et al., 2010).''') -parser.add_argument('-g','--gridspecdir', type=str, required=True, +parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') cmdLineArgs = parser.parse_args() rootGroup = netCDF4.Dataset( cmdLineArgs.monthly_file ) if 'MLD_003' not in rootGroup.variables: raise Exception('Could not find "MLD_003" in file "%s"'%(cmdLineArgs.monthly_file)) -x = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['x'][::2,::2] -y = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['y'][::2,::2] -msk = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_mask.nc').variables['mask'][:] -area = msk*netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) +if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') +if os.path.isdir(cmdLineArgs.gridspec): + x = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][::2,::2] + y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][::2,::2] + msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] + area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) +elif os.path.isfile(cmdLineArgs.gridspec): + x = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[::2,::2] + y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[::2,::2] + msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] + area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) +else: + raise ValueError('Unable to extract grid information from gridspec directory/tar file.') variable = rootGroup.variables['MLD_003'] shape = variable.shape diff --git a/tools/analysis/Makefile b/tools/analysis/Makefile new file mode 100644 index 0000000000..5a30e30caf --- /dev/null +++ b/tools/analysis/Makefile @@ -0,0 +1,47 @@ +# Makefile to test the analysis scripts called by frepp. + +# Settings for the 0.5 degree configuration +EXPDIR = /archive/John.Dunne/ulm_201505_awg_v20151106_mom6sis2_2015.12.17/CM4_c96L32_am4g7_2000_OMp5_DR_H3_MEKE/gfdl.ncrc2-intel-prod-openmp +GRIDSPEC = /archive/gold/datasets/OM4_05/mosaic.v20151203.unpacked +WOA = /archive/gold/datasets/OM4_05/obs/WOA05_ptemp_salt_annual.v2015.12.03.nc + +# Settings for the 0.25 degree configuration +#EXPDIR = /archive/Alistair.Adcroft/ulm_201505_awg_v20151106_mom6sis2_2016.01.16/CM4_c96L32_am4g7_2000_lowmix_hycom5_BS/gfdl.ncrc3-intel15-prod-openmp +#GRIDSPEC = /archive/gold/datasets/OM4_025/mosaic.v20140610.tar +#WOA = /archive/gold/datasets/OM4_025/obs/WOA05_ptemp_salt_annual.v20150310.nc + +all: sst sss zonalT zonalS zave sections moc mld heattransport + +sst: + ./SST_bias_WOA05.py -w $(WOA) -g $(GRIDSPEC) -l 0001-0005 $(EXPDIR)/pp/ocean_annual_z/av/annual_5yr/ocean_annual_z.0001-0005.ann.nc + +sss: + ./SSS_bias_WOA05.py -w $(WOA) -g $(GRIDSPEC) -l 0001-0005 $(EXPDIR)/pp/ocean_annual_z/av/annual_5yr/ocean_annual_z.0001-0005.ann.nc + +zonalT: + ./zonal_T_bias_WOA05.py -w $(WOA) -g $(GRIDSPEC) -l 0001-0005 $(EXPDIR)/pp/ocean_annual_z/av/annual_5yr/ocean_annual_z.0001-0005.ann.nc + +zonalS: + ./zonal_S_bias_WOA05.py -w $(WOA) -g $(GRIDSPEC) -l 0001-0005 $(EXPDIR)/pp/ocean_annual_z/av/annual_5yr/ocean_annual_z.0001-0005.ann.nc + +zave: + ./depth_average_T_bias.py -zb 100 -w $(WOA) -g $(GRIDSPEC) -l 0001-0005 $(EXPDIR)/pp/ocean_annual_z/av/annual_5yr/ocean_annual_z.0001-0005.ann.nc + ./depth_average_T_bias.py -zb 300 -w $(WOA) -g $(GRIDSPEC) -l 0001-0005 $(EXPDIR)/pp/ocean_annual_z/av/annual_5yr/ocean_annual_z.0001-0005.ann.nc + ./depth_average_T_bias.py -zb 700 -w $(WOA) -g $(GRIDSPEC) -l 0001-0005 $(EXPDIR)/pp/ocean_annual_z/av/annual_5yr/ocean_annual_z.0001-0005.ann.nc + ./depth_average_T_bias.py -zb 2000 -w $(WOA) -g $(GRIDSPEC) -l 0001-0005 $(EXPDIR)/pp/ocean_annual_z/av/annual_5yr/ocean_annual_z.0001-0005.ann.nc + ./depth_average_T_bias.py -zb 6500 -w $(WOA) -g $(GRIDSPEC) -l 0001-0005 $(EXPDIR)/pp/ocean_annual_z/av/annual_5yr/ocean_annual_z.0001-0005.ann.nc + +sections: + ./vertical_sections_annual_bias_WOA05.py -w $(WOA) -g $(GRIDSPEC) -l 0001-0005 $(EXPDIR)/pp/ocean_annual_z/av/annual_5yr/ocean_annual_z.0001-0005.ann.nc + +moc: + ./meridional_overturning.py -g $(GRIDSPEC) $(EXPDIR)/pp/ocean_annual_z/av/annual_5yr/ocean_annual_z.0001-0005.ann.nc + +mld: + ./MLD_003.py -g $(GRIDSPEC) $(EXPDIR)/pp/ocean_monthly/ts/monthly/5yr/ocean_monthly.000101-000512.MLD_003.nc + +heattransport: + ./poleward_heat_transport.py -l 0001-0005 -g $(GRIDSPEC) $(EXPDIR)/pp/ocean_annual/av/annual_5yr/ocean_annual.0001-0005.ann.nc + +clean: + rm -f *.png diff --git a/tools/analysis/SSS_bias_WOA05.py b/tools/analysis/SSS_bias_WOA05.py index 41db6c31db..3ef9cf6230 100755 --- a/tools/analysis/SSS_bias_WOA05.py +++ b/tools/analysis/SSS_bias_WOA05.py @@ -3,6 +3,8 @@ import netCDF4 import numpy import m6plot +import m6toolbox +import os import sys def run(): @@ -12,7 +14,7 @@ def run(): parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'salt'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') - parser.add_argument('-g','--gridspecdir', type=str, required=True, + parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') parser.add_argument('-w','--woa', type=str, required=True, help='''File containing WOA (or obs) data to compare against.''') @@ -21,12 +23,27 @@ def run(): def main(cmdLineArgs,stream=None): numpy.seterr(divide='ignore', invalid='ignore', over='ignore') # To avoid warnings + + if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') + if os.path.isdir(cmdLineArgs.gridspec): + x = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][::2,::2] + xcenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] + y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][::2,::2] + ycenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] + msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] + area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_topog.nc').variables['depth'][:] + elif os.path.isfile(cmdLineArgs.gridspec): + x = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[::2,::2] + xcenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[1::2,1::2] + y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[::2,::2] + ycenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2] + msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] + area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_topog.nc','depth')[:] + else: + raise ValueError('Unable to extract grid information from gridspec directory/tar file.') - x = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['x'][::2,::2] - y = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['y'][::2,::2] - msk = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_mask.nc').variables['mask'][:] - area = msk*netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) - msk = numpy.ma.array(msk, mask=(msk==0)) Sobs = netCDF4.Dataset( cmdLineArgs.woa ).variables['salt'] if len(Sobs.shape)==3: Sobs = Sobs[0] diff --git a/tools/analysis/SST_bias_WOA05.py b/tools/analysis/SST_bias_WOA05.py index 3d25b166e2..6a4b18ee04 100755 --- a/tools/analysis/SST_bias_WOA05.py +++ b/tools/analysis/SST_bias_WOA05.py @@ -3,6 +3,8 @@ import netCDF4 import numpy import m6plot +import m6toolbox +import os import sys def run(): @@ -12,7 +14,7 @@ def run(): parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'temp'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') - parser.add_argument('-g','--gridspecdir', type=str, required=True, + parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') parser.add_argument('-w','--woa', type=str, required=True, help='''File containing WOA (or obs) data to compare against.''') @@ -21,12 +23,26 @@ def run(): def main(cmdLineArgs,stream=None): numpy.seterr(divide='ignore', invalid='ignore', over='ignore') # To avoid warnings - - x = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['x'][::2,::2] - y = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['y'][::2,::2] - msk = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_mask.nc').variables['mask'][:] - area = msk*netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) - msk = numpy.ma.array(msk, mask=(msk==0)) + + if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') + if os.path.isdir(cmdLineArgs.gridspec): + x = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][::2,::2] + xcenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] + y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][::2,::2] + ycenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] + msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] + area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_topog.nc').variables['depth'][:] + elif os.path.isfile(cmdLineArgs.gridspec): + x = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[::2,::2] + xcenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[1::2,1::2] + y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[::2,::2] + ycenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2] + msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] + area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_topog.nc','depth')[:] + else: + raise ValueError('Unable to extract grid information from gridspec directory/tar file.') Tobs = netCDF4.Dataset( cmdLineArgs.woa ) if 'temp' in Tobs.variables: Tobs = Tobs.variables['temp'] diff --git a/tools/analysis/basin_test.py b/tools/analysis/basin_test.py new file mode 100644 index 0000000000..88c662bf03 --- /dev/null +++ b/tools/analysis/basin_test.py @@ -0,0 +1,54 @@ +import netCDF4 +import numpy +import m6plot +import m6toolbox +import matplotlib.pyplot as plt +import matplotlib.colors +import matplotlib.patches as mpatches +import os +import sys + +gridspec = '/archive/gold/datasets/OM4_025/mosaic.v20140610.unpacked' + +x = netCDF4.Dataset(gridspec+'/ocean_hgrid.nc').variables['x'][::2,::2] +xcenter = netCDF4.Dataset(gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] +y = netCDF4.Dataset(gridspec+'/ocean_hgrid.nc').variables['y'][::2,::2] +ycenter = netCDF4.Dataset(gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] +msk = netCDF4.Dataset(gridspec+'/ocean_mask.nc').variables['mask'][:] +area = msk*netCDF4.Dataset(gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) +depth = netCDF4.Dataset(gridspec+'/ocean_topog.nc').variables['depth'][:] + +code = m6toolbox.genBasinMasks(xcenter, ycenter, depth) + + +cmap = matplotlib.colors.ListedColormap(['#a1a1a1','#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99', + '#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a']) + +colors = ['#a1a1a1','#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99', + '#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a'] + +cmap = matplotlib.colors.ListedColormap(colors) + +plt.figure() +plt.pcolormesh(code,cmap=cmap) +plt.title('Basin Mask') + +ax = plt.gca() +proxy1 = plt.Rectangle((0, 0), 1, 1, fc=colors[1]) +proxy2 = plt.Rectangle((0, 0), 1, 1, fc=colors[2]) +proxy3 = plt.Rectangle((0, 0), 1, 1, fc=colors[3]) +proxy4 = plt.Rectangle((0, 0), 1, 1, fc=colors[4]) +proxy5 = plt.Rectangle((0, 0), 1, 1, fc=colors[5]) +proxy6 = plt.Rectangle((0, 0), 1, 1, fc=colors[6]) +proxy7 = plt.Rectangle((0, 0), 1, 1, fc=colors[7]) +proxy8 = plt.Rectangle((0, 0), 1, 1, fc=colors[8]) +proxy9 = plt.Rectangle((0, 0), 1, 1, fc=colors[9]) +proxy10 = plt.Rectangle((0, 0), 1, 1, fc=colors[10]) + +labels = ['Southern Ocean (1)','Atlantic Ocean (2)','Pacific Ocean (3)','Arctic Ocean (4)', + 'Indian Ocean (5)','Mediterranean Sea (6)','Black Sea (7)','Hudson Bay (8)', + 'Baltic Sea (9)','Red Sea (10)'] +proxies = [proxy1,proxy2,proxy3,proxy4,proxy5,proxy6,proxy7,proxy8,proxy9,proxy10] +ax.legend(proxies,labels,bbox_to_anchor=(0.01, -0.20), loc=2, borderaxespad=0., ncol=2) + + diff --git a/tools/analysis/depth_average_T_bias.py b/tools/analysis/depth_average_T_bias.py index 87b97d814d..397ba9ba65 100755 --- a/tools/analysis/depth_average_T_bias.py +++ b/tools/analysis/depth_average_T_bias.py @@ -1,7 +1,9 @@ #!/usr/bin/env python +import m6toolbox import netCDF4 import numpy +import os import m6plot try: import argparse @@ -11,7 +13,7 @@ parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'temp' and 'e'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') -parser.add_argument('-g','--gridspecdir', type=str, required=True, +parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') parser.add_argument('-w','--woa', type=str, required=True, help='''File containing WOA (or obs) data to compare against.''') @@ -19,11 +21,26 @@ parser.add_argument('-zb','--zbottom', type=float, default=300., help='''Depth (+ve) over which to end the average (default 300).''') cmdLineArgs = parser.parse_args() -x = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['x'][::2,::2] -y = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['y'][::2,::2] -msk = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_mask.nc').variables['mask'][:] -area = msk*netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) -depth = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_topog.nc').variables['depth'][:] +if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') +if os.path.isdir(cmdLineArgs.gridspec): + x = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][::2,::2] + xcenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] + y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][::2,::2] + ycenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] + msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] + area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_topog.nc').variables['depth'][:] +elif os.path.isfile(cmdLineArgs.gridspec): + x = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[::2,::2] + xcenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[1::2,1::2] + y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[::2,::2] + ycenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2] + msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] + area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_topog.nc','depth')[:] +else: + raise ValueError('Unable to extract grid information from gridspec directory/tar file.') + depth = numpy.ma.array(depth, mask=depth<=cmdLineArgs.ztop) lDepth = cmdLineArgs.zbottom uDepth = cmdLineArgs.ztop diff --git a/tools/analysis/m6toolbox.py b/tools/analysis/m6toolbox.py index 3427174d01..3643c92776 100644 --- a/tools/analysis/m6toolbox.py +++ b/tools/analysis/m6toolbox.py @@ -2,6 +2,8 @@ A collection of useful functions... """ import numpy as np +import tarfile +from scipy.io import netcdf def section2quadmesh(x, z, q, representation='pcm'): """ @@ -127,12 +129,15 @@ def ice9(i, j, source, xcyclic=True, tripolar=True): if i>0: stack.add( (j,i-1) ) elif xcyclic: stack.add( (j,ni-1) ) if i0: stack.add( (j-1,i) ) if j-zCellTop] = 1 return wet +def nearestJI(x, y, (x0, y0)): + """ + Find (j,i) of cell with center nearest to (x0,y0). + """ + return np.unravel_index( ((x-x0)**2 + (y-y0)**2).argmin() , x.shape) + +def readNCFromTar(tar,file,var): + TF = tarfile.open(tar,'r') + member = [m for m in TF.getmembers() if file in m.name][0] + nc = netcdf.netcdf_file(TF.extractfile(member),'r') + return nc.variables[var] + +def southOf(x, y, xy0, xy1): + """ + Returns 1 for point south/east of the line that passes through xy0-xy1, 0 otherwise. + """ + x0 = xy0[0]; y0 = xy0[1]; x1 = xy1[0]; y1 = xy1[1] + dx = x1 - x0; dy = y1 - y0 + Y = (x-x0)*dy - (y-y0)*dx + Y[Y>=0] = 1; Y[Y<=0] = 0 + return Y + +def genBasinMasks(x,y,depth,verbose=False): + if verbose: print 'Generating global wet mask ...', + wet = ice9Wrapper(x, y, depth, (0,-35)) # All ocean points seeded from South Atlantic + if verbose: print 'done.' + + code = 0*wet + + if verbose: print 'Finding Cape of Good Hope ...', + tmp = 1 - wet; tmp[x<-30] = 0 + tmp = ice9Wrapper(x, y, tmp, (20,-30.)) + yCGH = (tmp*y).min() + if verbose: print 'done.', yCGH + + if verbose: print 'Finding Melbourne ...', + tmp = 1 - wet; tmp[x>-180] = 0 + tmp = ice9Wrapper(x, y, tmp, (-220,-25.)) + yMel = (tmp*y).min() + if verbose: print 'done.', yMel + + if verbose: print 'Processing Persian Gulf ...' + tmp = wet*( 1-southOf(x, y, (55.,23.), (56.5,27.)) ) + tmp = ice9Wrapper(x, y, tmp, (53.,25.)) + code[tmp>0] = 11 + wet = wet - tmp # Removed named points + + if verbose: print 'Processing Red Sea ...' + tmp = wet*( 1-southOf(x, y, (40.,11.), (45.,13.)) ) + tmp = ice9Wrapper(x, y, tmp, (40.,18.)) + code[tmp>0] = 10 + wet = wet - tmp # Removed named points + + if verbose: print 'Processing Black Sea ...' + tmp = wet*( 1-southOf(x, y, (26.,42.), (32.,40.)) ) + tmp = ice9Wrapper(x, y, tmp, (32.,43.)) + code[tmp>0] = 7 + wet = wet - tmp # Removed named points + + if verbose: print 'Processing Mediterranean ...' + tmp = wet*( southOf(x, y, (-5.7,35.5), (-5.7,36.5)) ) + tmp = ice9Wrapper(x, y, tmp, (4.,38.)) + code[tmp>0] = 6 + wet = wet - tmp # Removed named points + + if verbose: print 'Processing Baltic ...' + tmp = wet*( southOf(x, y, (8.6,56.), (8.6,60.)) ) + tmp = ice9Wrapper(x, y, tmp, (10.,58.)) + code[tmp>0] = 9 + wet = wet - tmp # Removed named points + + if verbose: print 'Processing Hudson Bay ...' + tmp = wet*( + ( 1-(1-southOf(x, y, (-95.,66.), (-83.5,67.5))) + *(1-southOf(x, y, (-83.5,67.5), (-84.,71.))) + )*( 1-southOf(x, y, (-70.,58.), (-70.,65.)) ) ) + tmp = ice9Wrapper(x, y, tmp, (-85.,60.)) + code[tmp>0] = 8 + wet = wet - tmp # Removed named points + + if verbose: print 'Processing Arctic ...' + tmp = wet*( + (1-southOf(x, y, (-171.,66.), (-166.,65.5))) * (1-southOf(x, y, (-64.,66.4), (-50.,68.5))) # Lab Sea + + southOf(x, y, (-50.,0.), (-50.,90.)) * (1- southOf(x, y, (0.,65.5), (360.,65.5)) ) # Denmark Strait + + southOf(x, y, (-18.,0.), (-18.,65.)) * (1- southOf(x, y, (0.,64.9), (360.,64.9)) ) # Iceland-Sweden + + southOf(x, y, (20.,0.), (20.,90.)) # Barents Sea + + (1-southOf(x, y, (-280.,55.), (-200.,65.))) + ) + tmp = ice9Wrapper(x, y, tmp, (0.,85.)) + code[tmp>0] = 4 + wet = wet - tmp # Removed named points + + if verbose: print 'Processing Pacific ...' + tmp = wet*( (1-southOf(x, y, (0.,yMel), (360.,yMel))) + -southOf(x, y, (-257,1), (-257,0))*southOf(x, y, (0,3), (1,3)) + -southOf(x, y, (-254.25,1), (-254.25,0))*southOf(x, y, (0,-5), (1,-5)) + -southOf(x, y, (-243.7,1), (-243.7,0))*southOf(x, y, (0,-8.4), (1,-8.4)) + -southOf(x, y, (-234.5,1), (-234.5,0))*southOf(x, y, (0,-8.9), (1,-8.9)) + ) + tmp = ice9Wrapper(x, y, tmp, (-150.,0.)) + code[tmp>0] = 3 + wet = wet - tmp # Removed named points + + if verbose: print 'Processing Atlantic ...' + tmp = wet*(1-southOf(x, y, (0.,yCGH), (360.,yCGH))) + tmp = ice9Wrapper(x, y, tmp, (-20.,0.)) + code[tmp>0] = 2 + wet = wet - tmp # Removed named points + + if verbose: print 'Processing Indian ...' + tmp = wet*(1-southOf(x, y, (0.,yCGH), (360.,yCGH))) + tmp = ice9Wrapper(x, y, tmp, (55.,0.)) + code[tmp>0] = 5 + wet = wet - tmp # Removed named points + + if verbose: print 'Processing Southern Ocean ...' + tmp = ice9Wrapper(x, y, wet, (0.,-55.)) + code[tmp>0] = 1 + wet = wet - tmp # Removed named points + + code[wet>0] = -9 + (j,i) = np.unravel_index( wet.argmax(), x.shape) + if j: + if verbose: print 'There are leftover points unassigned to a basin code' + while j: + print x[j,i],y[j,i],[j,i] + wet[j,i]=0 + (j,i) = np.unravel_index( wet.argmax(), x.shape) + else: + if verbose: print 'All points assigned a basin code' + + if verbose: + print """ +Basin codes: +----------------------------------------------------------- + (1) Southern Ocean (6) Mediterranean Sea + (2) Atlantic Ocean (7) Black Sea + (3) Pacific Ocean (8) Hudson Bay + (4) Arctic Ocean (9) Baltic Sea + (5) Indian Ocean (10) Red Sea + + """ + + return code + + # Tests if __name__ == '__main__': @@ -182,3 +333,4 @@ def maskFromDepth(depth, zCellTop): plt.pcolormesh(X, Z, Q) plt.show() + diff --git a/tools/analysis/meridional_overturning.py b/tools/analysis/meridional_overturning.py index bcf53f012f..cc069df0da 100755 --- a/tools/analysis/meridional_overturning.py +++ b/tools/analysis/meridional_overturning.py @@ -3,7 +3,9 @@ import netCDF4 import numpy import m6plot +import m6toolbox import matplotlib.pyplot as plt +import os import sys def run(): @@ -13,19 +15,36 @@ def run(): parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'vh' and 'e'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') - parser.add_argument('-g','--gridspecdir', type=str, required=True, + parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') cmdLineArgs = parser.parse_args() main(cmdLineArgs) def main(cmdLineArgs,stream=None): - #x = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['x'][::2,::2] - y = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['y'][::2,::2] - msk = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_mask.nc').variables['mask'][:] - area = msk*netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) - depth = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_topog.nc').variables['depth'][:] - basin_code = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/basin_codes.v20140629.nc').variables['basin'][:] - + if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') + if os.path.isdir(cmdLineArgs.gridspec): + x = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][::2,::2] + xcenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] + y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][::2,::2] + ycenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] + msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] + area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_topog.nc').variables['depth'][:] + try: basin_code = netCDF4.Dataset(cmdLineArgs.gridspec+'/basin_codes.nc').variables['basin'][:] + except: basin_code = m6toolbox.genBasinMasks(xcenter, ycenter, depth) + elif os.path.isfile(cmdLineArgs.gridspec): + x = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[::2,::2] + xcenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[1::2,1::2] + y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[::2,::2] + ycenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2] + msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] + area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_topog.nc','depth')[:] + try: basin_code = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'basin_codes.nc','basin')[:] + except: basin_code = m6toolbox.genBasinMasks(xcenter, ycenter, depth) + else: + raise ValueError('Unable to extract grid information from gridspec directory/tar file.') + if stream != None: if len(stream) != 2: raise ValueError('If specifying output streams, exactly two streams are needed for this analysis') @@ -36,8 +55,10 @@ def main(cmdLineArgs,stream=None): elif 'vmo' in rootGroup.variables: varName = 'vmo'; conversion_factor = 1.e-9 else: raise Exception('Could not find "vh" or "vmo" in file "%s"'%(cmdLineArgs.annual_file)) - if len(rootGroup.variables[varName].shape)==4: VHmod = rootGroup.variables[varName][:].mean(axis=0).filled(0.) - else: VHmod = rootGroup.variables[varName][:].filled(0.) + if len(rootGroup.variables[varName].shape)==4: VHmod = rootGroup.variables[varName][:].mean(axis=0) + else: VHmod = rootGroup.variables[varName][:] + try: VHmod = VHmod.filled(0.) + except: pass if 'e' in rootGroup.variables: Zmod = rootGroup.variables['e'][0] elif 'zw' in rootGroup.variables: zw = rootGroup.variables['zw'][:] @@ -99,7 +120,7 @@ def findExtrema(y, z, psi, min_lat=-90., max_lat=90., min_depth=0., mult=1.): plt.clf() m = 0*basin_code; m[(basin_code==2) | (basin_code==4) | (basin_code==6) | (basin_code==7) | (basin_code==8)]=1 ci=m6plot.pmCI(0.,22.,2.) - z = (m*Zmod).min(axis=-1); psiPlot = MOCpsi(VHmod, vmsk=m*numpy.roll(m,1,axis=1))*conversion_factor + z = (m*Zmod).min(axis=-1); psiPlot = MOCpsi(VHmod, vmsk=m*numpy.roll(m,-1,axis=-2))*conversion_factor yy = y[1:,:].max(axis=-1)+0*z plotPsi(yy, z, psiPlot, ci, 'Atlantic MOC [Sv]') plt.xlabel(r'Latitude [$\degree$N]') diff --git a/tools/analysis/poleward_heat_transport.py b/tools/analysis/poleward_heat_transport.py index 8a05209921..ccbcb5df64 100755 --- a/tools/analysis/poleward_heat_transport.py +++ b/tools/analysis/poleward_heat_transport.py @@ -3,7 +3,9 @@ import netCDF4 import numpy import m6plot +import m6toolbox import matplotlib.pyplot as plt +import os import sys import warnings @@ -14,22 +16,39 @@ def run(): parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'T_ady_2d' and 'T_diffy_2d'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') - parser.add_argument('-g','--gridspecdir', type=str, required=True, - help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') + parser.add_argument('-g','--gridspec', type=str, required=True, + help='''Directory or tarfile containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') cmdLineArgs = parser.parse_args() main(cmdLineArgs) def main(cmdLineArgs,stream=None): - #x = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['x'][::2,::2] - y = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['y'][::2,::2] - msk = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_mask.nc').variables['mask'][:] - area = msk*netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) - depth = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_topog.nc').variables['depth'][:] - basin_code = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/basin_codes.v20140629.nc').variables['basin'][:] - + if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') + if os.path.isdir(cmdLineArgs.gridspec): + x = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][::2,::2] + xcenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] + y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][::2,::2] + ycenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] + msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] + area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_topog.nc').variables['depth'][:] + try: basin_code = netCDF4.Dataset(cmdLineArgs.gridspec+'/basin_codes.nc').variables['basin'][:] + except: basin_code = m6toolbox.genBasinMasks(xcenter, ycenter, depth) + elif os.path.isfile(cmdLineArgs.gridspec): + x = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[::2,::2] + xcenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[1::2,1::2] + y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[::2,::2] + ycenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2] + msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] + area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_topog.nc','depth')[:] + try: basin_code = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'basin_codes.nc','basin')[:] + except: basin_code = m6toolbox.genBasinMasks(xcenter, ycenter, depth) + else: + raise ValueError('Unable to extract grid information from gridspec directory/tar file.') + if stream != None: - if len(stream) != 2: - raise ValueError('If specifying output streams, exactly two streams are needed for this analysis') + if len(stream) != 3: + raise ValueError('If specifying output streams, exactly three streams are needed for this analysis') rootGroup = netCDF4.MFDataset( cmdLineArgs.annual_file ) if 'T_ady_2d' in rootGroup.variables: @@ -79,12 +98,11 @@ def annotatePlot(label): plt.savefig(stream[0]) else: plt.savefig(cmdLineArgs.outdir+'/HeatTransport_global.png') - plt.show(block=False) # Atlantic Heat Transport plt.clf() - m = 0*basin_code; m[(basin_code==2) | (basin_code==4)] = 1 - HTplot = heatTrans(advective, diffusive, vmask=m*numpy.roll(m,1,axis=1)) + m = 0*basin_code; m[(basin_code==2) | (basin_code==4) | (basin_code==6) | (basin_code==7) | (basin_code==8)] = 1 + HTplot = heatTrans(advective, diffusive, vmask=m*numpy.roll(m,-1,axis=-2)) yy = y[1:,:].max(axis=-1) HTplot[yy<-34] = numpy.nan plotHeatTrans(yy,HTplot,title='Atlantic Y-Direction Heat Transport [PW]') @@ -92,15 +110,14 @@ def annotatePlot(label): plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) if diffusive == None: annotatePlot('Warning: Diffusive component of transport is missing.') if stream != None: - plt.savefig(stream[0]) + plt.savefig(stream[1]) else: plt.savefig(cmdLineArgs.outdir+'/HeatTransport_Atlantic.png') - plt.show(block=False) # Indo-Pacific Heat Transport plt.clf() m = 0*basin_code; m[(basin_code==3) | (basin_code==5)] = 1 - HTplot = heatTrans(advective, diffusive, vmask=m*numpy.roll(m,1,axis=1)) + HTplot = heatTrans(advective, diffusive, vmask=m*numpy.roll(m,-1,axis=-2)) yy = y[1:,:].max(axis=-1) HTplot[yy<-38] = numpy.nan plotHeatTrans(yy,HTplot,title='Indo-Pacific Y-Direction Heat Transport [PW]') @@ -108,10 +125,9 @@ def annotatePlot(label): if diffusive == None: annotatePlot('Warning: Diffusive component of transport is missing.') plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) if stream != None: - plt.savefig(stream[0]) + plt.savefig(stream[2]) else: plt.savefig(cmdLineArgs.outdir+'/HeatTransport_IndoPac.png') - plt.show(block=False) if __name__ == '__main__': run() diff --git a/tools/analysis/vertical_sections_annual_bias_WOA05.py b/tools/analysis/vertical_sections_annual_bias_WOA05.py index 37b042127b..dc4f6b56f9 100755 --- a/tools/analysis/vertical_sections_annual_bias_WOA05.py +++ b/tools/analysis/vertical_sections_annual_bias_WOA05.py @@ -3,6 +3,8 @@ import netCDF4 import numpy import m6plot +import m6toolbox +import os import matplotlib.pyplot as plt try: import argparse @@ -12,19 +14,34 @@ parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'temp' and 'e'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') -parser.add_argument('-g','--gridspecdir', type=str, required=True, +parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') parser.add_argument('-w','--woa', type=str, required=True, help='''File containing WOA (or obs) data to compare against.''') cmdLineArgs = parser.parse_args() -x = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['x'][1::2,1::2] -y = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['y'][1::2,1::2] -xg = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['x'][0::2,0::2] -yg = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['y'][0::2,0::2] -msk = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_mask.nc').variables['mask'][:] -area = msk*netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) -basin = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/basin_codes.nc').variables['basin'][:] +if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') +if os.path.isdir(cmdLineArgs.gridspec): + x = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] + xg = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][0::2,0::2] + y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] + yg = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][0::2,0::2] + msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] + area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_topog.nc').variables['depth'][:] +elif os.path.isfile(cmdLineArgs.gridspec): + x = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[1::2,1::2] + xg = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[0::2,0::2] + y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2] + yg = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[0::2,0::2] + msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] + area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_topog.nc','depth')[:] +else: + raise ValueError('Unable to extract grid information from gridspec directory/tar file.') + +# Basin codes +basin = m6toolbox.genBasinMasks(x, y, depth) # All ocean points seeded from South Atlantic obsRootGroup = netCDF4.Dataset( cmdLineArgs.woa ) if 'temp' in obsRootGroup.variables: OTvar = 'temp' diff --git a/tools/analysis/zonal_S_bias_WOA05.py b/tools/analysis/zonal_S_bias_WOA05.py index f4573384b1..a43c718576 100755 --- a/tools/analysis/zonal_S_bias_WOA05.py +++ b/tools/analysis/zonal_S_bias_WOA05.py @@ -1,6 +1,8 @@ #!/usr/bin/env python import netCDF4 +import os +import m6toolbox import numpy import m6plot @@ -11,16 +13,36 @@ parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'salt' and 'e'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') -parser.add_argument('-g','--gridspecdir', type=str, required=True, +parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') parser.add_argument('-w','--woa', type=str, required=True, help='''File containing WOA (or obs) data to compare against.''') cmdLineArgs = parser.parse_args() -y = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['y'][1::2,1::2].max(axis=-1) -msk = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_mask.nc').variables['mask'][:] -area = msk*netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) -basin = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/basin_codes.nc').variables['basin'][:] +if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') +if os.path.isdir(cmdLineArgs.gridspec): + xcenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] + y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2].max(axis=-1) + ycenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] + msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] + area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_topog.nc').variables['depth'][:] + try: basin = netCDF4.Dataset(cmdLineArgs.gridspec+'/basin_codes.nc').variables['basin'][:] + except: basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) +elif os.path.isfile(cmdLineArgs.gridspec): + xcenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[1::2,1::2] + y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2].max(axis=-1) + ycenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2] + msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] + area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_topog.nc','depth')[:] + try: basin = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'basin_codes.nc','basin')[:] + except: basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) +else: + raise ValueError('Unable to extract grid information from gridspec directory/tar file.') + +# Basin codes +basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) # All ocean points seeded from South Atlantic Sobs = netCDF4.Dataset( cmdLineArgs.woa ).variables['salt'] if len(Sobs.shape)==3: Sobs = Sobs[:] diff --git a/tools/analysis/zonal_T_bias_WOA05.py b/tools/analysis/zonal_T_bias_WOA05.py index 5cbebf6ecb..6b7b36fc3f 100755 --- a/tools/analysis/zonal_T_bias_WOA05.py +++ b/tools/analysis/zonal_T_bias_WOA05.py @@ -2,6 +2,8 @@ import netCDF4 import numpy +import m6toolbox +import os import m6plot try: import argparse @@ -11,16 +13,33 @@ parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'temp' and 'e'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') -parser.add_argument('-g','--gridspecdir', type=str, required=True, +parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') parser.add_argument('-w','--woa', type=str, required=True, help='''File containing WOA (or obs) data to compare against.''') cmdLineArgs = parser.parse_args() -y = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['y'][1::2,1::2].max(axis=-1) -msk = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_mask.nc').variables['mask'][:] -area = msk*netCDF4.Dataset(cmdLineArgs.gridspecdir+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) -basin = netCDF4.Dataset(cmdLineArgs.gridspecdir+'/basin_codes.nc').variables['basin'][:] +if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') +if os.path.isdir(cmdLineArgs.gridspec): + xcenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] + y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2].max(axis=-1) + ycenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] + msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] + area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_topog.nc').variables['depth'][:] + try: basin = netCDF4.Dataset(cmdLineArgs.gridspec+'/basin_codes.nc').variables['basin'][:] + except: basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) +elif os.path.isfile(cmdLineArgs.gridspec): + xcenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[1::2,1::2] + y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2].max(axis=-1) + ycenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2] + msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] + area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_topog.nc','depth')[:] + try: basin = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'basin_codes.nc','basin')[:] + except: basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) +else: + raise ValueError('Unable to extract grid information from gridspec directory/tar file.') Tobs = netCDF4.Dataset( cmdLineArgs.woa ) if 'temp' in Tobs.variables: Tobs = Tobs.variables['temp'] From 12cabc9610cf6a471a544e13aba2be6ff42e7ea6 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Thu, 18 Feb 2016 13:49:16 -0500 Subject: [PATCH 3/9] Updated .frepp files - Added heat transport diagnostic to OM4_025 example - Added .frepp files for OM4_05 example --- ice_ocean_SIS2/OM4_025/ocn_annual.frepp | 3 + ice_ocean_SIS2/OM4_05/ocn_annual.frepp | 118 +++++++++++++++++++++ ice_ocean_SIS2/OM4_05/ocn_annual_z.frepp | 125 +++++++++++++++++++++++ ice_ocean_SIS2/OM4_05/ocn_monthly.frepp | 106 +++++++++++++++++++ 4 files changed, 352 insertions(+) create mode 100644 ice_ocean_SIS2/OM4_05/ocn_annual.frepp create mode 100644 ice_ocean_SIS2/OM4_05/ocn_annual_z.frepp create mode 100644 ice_ocean_SIS2/OM4_05/ocn_monthly.frepp diff --git a/ice_ocean_SIS2/OM4_025/ocn_annual.frepp b/ice_ocean_SIS2/OM4_025/ocn_annual.frepp index 1f3edf0cab..58c330c1a6 100644 --- a/ice_ocean_SIS2/OM4_025/ocn_annual.frepp +++ b/ice_ocean_SIS2/OM4_025/ocn_annual.frepp @@ -110,6 +110,9 @@ set mosaicDir=/archive/gold/datasets/OM4_025/mosaic.v20140610.unpacked echo The following command is going to be invoked +mkdir -p $out_dir/ocean_${yr1}-${yr2}/heat_transport +$script_dir/poleward_heat_transport.py -g $mosaicDir -o $out_dir/ocean_${yr1}-${yr2}/heat_transport -l ${yr1}-${yr2} $in_data_dir/ocean_annual.$yr1-$yr2.ann.nc + mkdir -p $out_dir/ocean_${yr1}-${yr2}/heat_salt_0_300m $script_dir/TS_depth_integrals.py -r $woa05 -s 0 -e 300 -g $mosaicDir -o $out_dir/ocean_${yr1}-${yr2}/heat_salt_0_300m -l ${yr1}-${yr2} $in_data_dir/ocean_annual.$yr1-$yr2.ann.nc \rm -f $out_dir/ocean_${yr1}-${yr2}/heat_salt_0_300m/*_heat_salt_0_300.nc # Clean up large nc files diff --git a/ice_ocean_SIS2/OM4_05/ocn_annual.frepp b/ice_ocean_SIS2/OM4_05/ocn_annual.frepp new file mode 100644 index 0000000000..8a98ce15f7 --- /dev/null +++ b/ice_ocean_SIS2/OM4_05/ocn_annual.frepp @@ -0,0 +1,118 @@ +#!/bin/csh -fx +#------------------------------------ +#PBS -N ocean_ts_annual +#PBS -l size=1 +#PBS -l walltime=04:00:00 +#PBS -r y +#PBS -j oe +#PBS -o +#PBS -q batch +#---------------------------------- + +echo $0 $* + +# ============== VARIABLES SET BY FREPP ============= +# original arguments passed to this script when it was created by frepp +set argu +# experiment name (as appears in output directory names) +set descriptor +# path to NetCDF files postprocessed by frepp, as specified in XML +set in_data_dir +# input data file[s], without any path prefix; +# currently only used by timeAverage diagnostics +set in_data_file +# final output directory for diagnostics generated by this script +set out_dir +# working directory -- do whatever you want in here +# we may have to create this (use "mkdir -p" to be sure) +# and then clean up at end +set WORKDIR +# a string to indicate the mode: "batch" or "interactive" +set mode +# actual start/end years of diagnostics (start_year & end_year in XML) +set yr1 +set yr2 +# alternate way to specify a single year (instead of yr1==yr2) +set specify_year +# Data years, only used in scripts using time series as input, to +# generate a Ferret descriptor file from consecutive NetCDF chunks. +# start year of first chunk +set databegyr +# end year of last chunk +set dataendyr +# chunk length (an integer number), as specified in XML +set datachunk +# atmospheric land mask file +set staticfile +set script_path +set fremodule +# a string: "monthly" or "annual" for timeseries data +set freq +# first year of integration (4-digits, e.g. the year of initial condition) +set MODEL_start_yr +# Specify MOM version; "om2" or "om3" because some files depend on mom's grid +set mom_version +# full path to the grid specification file, which contains the land/sea mask +set gridspecfile +# history directory where the original "*.nc.cpio" files are found +set hist_dir +set nlon +set nlat +set frexml +set freanalysismodule +set stdoutdir +set analysis_options +set platform +set target +# ============== END OF VARIABLES SET BY FREPP ============= +# Invoke script with arguments tcsh ocn_annual.frepp in_data_dir yr1 yr2 out_dir +if ($#argv == "4") then + if ($in_data_dir == "") set in_data_dir = $1 + if ($yr1 == "") set yr1 = $2 + if ($yr2 == "") set yr2 = $3 + if ($out_dir == "") set out_dir = $4 +endif + +#set freanalysismodule = fre-analysis/test + +# make sure valid platform and required modules are loaded +if (`gfdl_platform` == "hpcs-csc") then + source $MODULESHOME/init/csh + module use -a /home/fms/local/modulefiles #/usr/local/paida/Modules + module purge +# module load $fremodule +# module load $freanalysismodule +# module load gcc + module load netcdf/4.2 + module load python/2.7.3 + module load intel_compilers + module load git +else + echo "ERROR: invalid platform" + #exit 1 +endif + +# check again? +#if (! $?FRE_ANALYSIS_HOME) then +# echo "ERROR: environment variable FRE_ANALYSIS_HOME not set." +# exit 1 +#endif + +# Build MIDAS libs and set PYTHONPATH +(cd ${out_dir}/mom6/tools/python ; make frepp_local ) +setenv PYTHONPATH ${out_dir}/mom6/tools/python/local/lib/python + +# Run script +set src_dir = ${out_dir}/mom6/tools/analysis +set script_dir = ${out_dir}/mom6/tools/analysis +set woa05=/archive/gold/datasets/OM4_05/obs/WOA05_ptemp_salt_annual.v2015.12.03.nc +set mosaicDir=/archive/gold/datasets/OM4_05/mosaic.v20151203.unpacked + +echo The following command is going to be invoked + +mkdir -p $out_dir/ocean_${yr1}-${yr2}/heat_transport +$script_dir/poleward_heat_transport.py -g $mosaicDir -o $out_dir/ocean_${yr1}-${yr2}/heat_transport -l ${yr1}-${yr2} $in_data_dir/ocean_annual.$yr1-$yr2.ann.nc + +mkdir -p $out_dir/ocean_${yr1}-${yr2}/heat_salt_0_300m +$script_dir/TS_depth_integrals.py -r $woa05 -s 0 -e 300 -g $mosaicDir -o $out_dir/ocean_${yr1}-${yr2}/heat_salt_0_300m -l ${yr1}-${yr2} $in_data_dir/ocean_annual.$yr1-$yr2.ann.nc +\rm -f $out_dir/ocean_${yr1}-${yr2}/heat_salt_0_300m/*_heat_salt_0_300.nc # Clean up large nc files diff --git a/ice_ocean_SIS2/OM4_05/ocn_annual_z.frepp b/ice_ocean_SIS2/OM4_05/ocn_annual_z.frepp new file mode 100644 index 0000000000..6181ce653c --- /dev/null +++ b/ice_ocean_SIS2/OM4_05/ocn_annual_z.frepp @@ -0,0 +1,125 @@ +#!/bin/csh -fx +#------------------------------------ +#PBS -N ocean_ts_annual +#PBS -l size=1 +#PBS -l walltime=04:00:00 +#PBS -r y +#PBS -j oe +#PBS -o +#PBS -q batch +#---------------------------------- + +echo $0 $* + +# ============== VARIABLES SET BY FREPP ============= +# original arguments passed to this script when it was created by frepp +set argu +# experiment name (as appears in output directory names) +set descriptor +# path to NetCDF files postprocessed by frepp, as specified in XML +set in_data_dir +# input data file[s], without any path prefix; +# currently only used by timeAverage diagnostics +set in_data_file +# final output directory for diagnostics generated by this script +set out_dir +# working directory -- do whatever you want in here +# we may have to create this (use "mkdir -p" to be sure) +# and then clean up at end +set WORKDIR +# a string to indicate the mode: "batch" or "interactive" +set mode +# actual start/end years of diagnostics (start_year & end_year in XML) +set yr1 +set yr2 +# alternate way to specify a single year (instead of yr1==yr2) +set specify_year +# Data years, only used in scripts using time series as input, to +# generate a Ferret descriptor file from consecutive NetCDF chunks. +# start year of first chunk +set databegyr +# end year of last chunk +set dataendyr +# chunk length (an integer number), as specified in XML +set datachunk +# atmospheric land mask file +set staticfile +set script_path +set fremodule +# a string: "monthly" or "annual" for timeseries data +set freq +# first year of integration (4-digits, e.g. the year of initial condition) +set MODEL_start_yr +# Specify MOM version; "om2" or "om3" because some files depend on mom's grid +set mom_version +# full path to the grid specification file, which contains the land/sea mask +set gridspecfile +# history directory where the original "*.nc.cpio" files are found +set hist_dir +set nlon +set nlat +set frexml +set freanalysismodule +set stdoutdir +set analysis_options +set platform +set target +# ============== END OF VARIABLES SET BY FREPP ============= +# Invoke script with arguments tcsh ocn_annual_z.frepp in_data_dir yr1 yr2 out_dir +if ($#argv == "4") then + if ($in_data_dir == "") set in_data_dir = $1 + if ($yr1 == "") set yr1 = $2 + if ($yr2 == "") set yr2 = $3 + if ($out_dir == "") set out_dir = $4 +endif + +#set freanalysismodule = fre-analysis/test + +# make sure valid platform and required modules are loaded +if (`gfdl_platform` == "hpcs-csc") then + source $MODULESHOME/init/csh + module use -a /home/fms/local/modulefiles + module purge + module load $fremodule + module load $freanalysismodule + module load gcc + module load netcdf/4.2 + module load python/2.7.3 +else + echo "ERROR: invalid platform" + #exit 1 +endif + +# check again? +if (! $?FRE_ANALYSIS_HOME) then + echo "ERROR: environment variable FRE_ANALYSIS_HOME not set." + #exit 1 +endif + +# Run script +set src_dir = ${out_dir}/mom6/tools/analysis +set script_dir = ${out_dir}/mom6/tools/analysis +set woa05=/archive/gold/datasets/OM4_05/obs/WOA05_ptemp_salt_annual.v2015.12.03.nc +set mosaicDir=/archive/gold/datasets/OM4_05/mosaic.v20151203.unpacked + +echo The following command is going to be invoked + +mkdir -p $out_dir/ocean_${yr1}-${yr2}/ptemp +$script_dir/SST_bias_WOA05.py -w $woa05 -g $mosaicDir -o $out_dir/ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} $in_data_dir/ocean_annual_z.$yr1-$yr2.ann.nc +$script_dir/zonal_T_bias_WOA05.py -w $woa05 -g $mosaicDir -o $out_dir/ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} $in_data_dir/ocean_annual_z.$yr1-$yr2.ann.nc +$script_dir/depth_average_T_bias.py -w $woa05 -g $mosaicDir -zb 100 -o $out_dir/ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} $in_data_dir/ocean_annual_z.$yr1-$yr2.ann.nc +$script_dir/depth_average_T_bias.py -w $woa05 -g $mosaicDir -zb 300 -o $out_dir/ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} $in_data_dir/ocean_annual_z.$yr1-$yr2.ann.nc +$script_dir/depth_average_T_bias.py -w $woa05 -g $mosaicDir -zb 700 -o $out_dir/ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} $in_data_dir/ocean_annual_z.$yr1-$yr2.ann.nc +$script_dir/depth_average_T_bias.py -w $woa05 -g $mosaicDir -zb 2000 -o $out_dir/ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} $in_data_dir/ocean_annual_z.$yr1-$yr2.ann.nc +$script_dir/depth_average_T_bias.py -w $woa05 -g $mosaicDir -zt 2000 -zb 4000 -o $out_dir/ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} $in_data_dir/ocean_annual_z.$yr1-$yr2.ann.nc +$script_dir/depth_average_T_bias.py -w $woa05 -g $mosaicDir -zb 6500 -o $out_dir/ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} $in_data_dir/ocean_annual_z.$yr1-$yr2.ann.nc + +mkdir -p $out_dir/ocean_${yr1}-${yr2}/salinity +$script_dir/SSS_bias_WOA05.py -w $woa05 -g $mosaicDir -o $out_dir/ocean_${yr1}-${yr2}/salinity -l ${yr1}-${yr2} $in_data_dir/ocean_annual_z.$yr1-$yr2.ann.nc +$script_dir/zonal_S_bias_WOA05.py -w $woa05 -g $mosaicDir -o $out_dir/ocean_${yr1}-${yr2}/salinity -l ${yr1}-${yr2} $in_data_dir/ocean_annual_z.$yr1-$yr2.ann.nc + +mkdir -p $out_dir/ocean_${yr1}-${yr2}/sections +$script_dir/vertical_sections_annual_bias_WOA05.py -g $mosaicDir -o $out_dir/ocean_${yr1}-${yr2}/sections -w $woa05 -l ${yr1}-${yr2} $in_data_dir/ocean_annual_z.$yr1-$yr2.ann.nc + +mkdir -p $out_dir/ocean_${yr1}-${yr2}/MOC +$script_dir/meridional_overturning.py -g $mosaicDir -o $out_dir/ocean_${yr1}-${yr2}/MOC -l ${yr1}-${yr2} $in_data_dir/ocean_annual_z.$yr1-$yr2.ann.nc diff --git a/ice_ocean_SIS2/OM4_05/ocn_monthly.frepp b/ice_ocean_SIS2/OM4_05/ocn_monthly.frepp new file mode 100644 index 0000000000..ebf0666c7a --- /dev/null +++ b/ice_ocean_SIS2/OM4_05/ocn_monthly.frepp @@ -0,0 +1,106 @@ +#!/bin/csh -fx +#------------------------------------ +#PBS -N ocean_ts_annual +#PBS -l size=1 +#PBS -l walltime=04:00:00 +#PBS -r y +#PBS -j oe +#PBS -o +#PBS -q batch +#---------------------------------- + +echo $0 $* + +# ============== VARIABLES SET BY FREPP ============= +# original arguments passed to this script when it was created by frepp +set argu +# experiment name (as appears in output directory names) +set descriptor +# path to NetCDF files postprocessed by frepp, as specified in XML +set in_data_dir +# input data file[s], without any path prefix; +# currently only used by timeAverage diagnostics +set in_data_file +# final output directory for diagnostics generated by this script +set out_dir +# working directory -- do whatever you want in here +# we may have to create this (use "mkdir -p" to be sure) +# and then clean up at end +set WORKDIR +# a string to indicate the mode: "batch" or "interactive" +set mode +# actual start/end years of diagnostics (start_year & end_year in XML) +set yr1 +set yr2 +# alternate way to specify a single year (instead of yr1==yr2) +set specify_year +# Data years, only used in scripts using time series as input, to +# generate a Ferret descriptor file from consecutive NetCDF chunks. +# start year of first chunk +set databegyr +# end year of last chunk +set dataendyr +# chunk length (an integer number), as specified in XML +set datachunk +# atmospheric land mask file +set staticfile +set script_path +set fremodule +# a string: "monthly" or "annual" for timeseries data +set freq +# first year of integration (4-digits, e.g. the year of initial condition) +set MODEL_start_yr +# Specify MOM version; "om2" or "om3" because some files depend on mom's grid +set mom_version +# full path to the grid specification file, which contains the land/sea mask +set gridspecfile +# history directory where the original "*.nc.cpio" files are found +set hist_dir +set nlon +set nlat +set frexml +set freanalysismodule +set stdoutdir +set analysis_options +set platform +set target +# ============== END OF VARIABLES SET BY FREPP ============= +# Invoke script with arguments tcsh ocn_monthly.frepp in_data_dir yr1 yr2 out_dir +if ($#argv == "4") then + if ($in_data_dir == "") set in_data_dir = $1 + if ($yr1 == "") set yr1 = $2 + if ($yr2 == "") set yr2 = $3 + if ($out_dir == "") set out_dir = $4 +endif + +#set freanalysismodule = fre-analysis/test + +# make sure valid platform and required modules are loaded +if (`gfdl_platform` == "hpcs-csc") then + source $MODULESHOME/init/csh + module use -a /home/fms/local/modulefiles + module purge + module load $fremodule + module load $freanalysismodule + module load gcc + module load netcdf/4.2 + module load python/2.7.3 +else + echo "ERROR: invalid platform" + #exit 1 +endif + +# check again? +if (! $?FRE_ANALYSIS_HOME) then + echo "ERROR: environment variable FRE_ANALYSIS_HOME not set." + #exit 1 +endif + +# Run script +set script_dir = ${out_dir}/mom6/tools/analysis +set mosaicDir=/archive/gold/datasets/OM4_05/mosaic.v20151203.unpacked + +echo The following command is going to be invoked + +mkdir -p $out_dir/ocean_${yr1}-${yr2}/Hosoda_MLD +$script_dir/MLD_003.py -g $mosaicDir -o $out_dir/ocean_${yr1}-${yr2}/Hosoda_MLD -l ${yr1}-${yr2} $in_data_dir/ocean_monthly.${yr1}01-${yr2}12.MLD_003.nc From ffe97247db9f843cf9f94675a79c27729c7e927f Mon Sep 17 00:00:00 2001 From: John Krasting Date: Thu, 18 Feb 2016 15:27:12 -0500 Subject: [PATCH 4/9] Bring ESM4 SIS2 up-to-date --- src/SIS2 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SIS2 b/src/SIS2 index a028d6679b..cb45714190 160000 --- a/src/SIS2 +++ b/src/SIS2 @@ -1 +1 @@ -Subproject commit a028d6679b7debb2e50b9e5952e4482f97177031 +Subproject commit cb457141901bd6367a97e9ed62307c254285991b From ea37d37cbd5556db544e180135d4a7b7ac02a35e Mon Sep 17 00:00:00 2001 From: John Krasting Date: Thu, 18 Feb 2016 15:52:34 -0500 Subject: [PATCH 5/9] Generic tracer update from JGJ - See MOM6 commit ae530a2e --- src/MOM6 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MOM6 b/src/MOM6 index e68537d2ad..ae530a2ef2 160000 --- a/src/MOM6 +++ b/src/MOM6 @@ -1 +1 @@ -Subproject commit e68537d2ad674c573e3b857d7afa8776b43034d2 +Subproject commit ae530a2ef2031b188351a0503b7285fb24f19cfe From f9b8b0ce790388fc2cd207bfa9c601d2952f6aef Mon Sep 17 00:00:00 2001 From: John Krasting Date: Tue, 1 Mar 2016 09:03:13 -0500 Subject: [PATCH 6/9] Updates for dora web analysis - Functionality for stringIO buffer streams - Added ocean heat transport observations --- tools/analysis/poleward_heat_transport.py | 76 ++++++- tools/analysis/zonal_S_bias_WOA05.py | 261 +++++++++++---------- tools/analysis/zonal_T_bias_WOA05.py | 264 ++++++++++++---------- 3 files changed, 352 insertions(+), 249 deletions(-) diff --git a/tools/analysis/poleward_heat_transport.py b/tools/analysis/poleward_heat_transport.py index ccbcb5df64..9bb97ddc38 100755 --- a/tools/analysis/poleward_heat_transport.py +++ b/tools/analysis/poleward_heat_transport.py @@ -9,6 +9,7 @@ import sys import warnings + def run(): try: import argparse except: raise Exception('This version of python is not new enough. python 2.7 or newer is required.') @@ -53,12 +54,12 @@ def main(cmdLineArgs,stream=None): rootGroup = netCDF4.MFDataset( cmdLineArgs.annual_file ) if 'T_ady_2d' in rootGroup.variables: varName = 'T_ady_2d' - if len(rootGroup.variables[varName].shape)==4: advective = rootGroup.variables[varName][:].mean(axis=0).filled(0.) + if len(rootGroup.variables[varName].shape)==3: advective = rootGroup.variables[varName][:].mean(axis=0).filled(0.) else: advective = rootGroup.variables[varName][:].filled(0.) else: raise Exception('Could not find "T_ady_2d" in file "%s"'%(cmdLineArgs.annual_file)) if 'T_diffy_2d' in rootGroup.variables: varName = 'T_diffy_2d' - if len(rootGroup.variables[varName].shape)==4: diffusive = rootGroup.variables[varName][:].mean(axis=0).filled(0.) + if len(rootGroup.variables[varName].shape)==3: diffusive = rootGroup.variables[varName][:].mean(axis=0).filled(0.) else: diffusive = rootGroup.variables[varName][:].filled(0.) else: diffusive = None @@ -74,25 +75,76 @@ def heatTrans(advective, diffusive=None, vmask=None): HT = HT.sum(axis=-1); HT = HT.squeeze() # sum in x-direction return HT - def plotHeatTrans(y, HT, title): + def plotHeatTrans(y, HT, title, xlim=(-80,90)): plt.plot(y, y*0., 'k', linewidth=0.5) - plt.plot(y, HT, 'k', linewidth=1.5) - plt.xlim(-80,90); plt.ylim(-2.0,3.0) + plt.plot(y, HT, 'r', linewidth=1.5,label='Model') + plt.xlim(xlim); plt.ylim(-2.5,3.0) plt.title(title) plt.grid(True) def annotatePlot(label): fig = plt.gcf() - fig.text(0.1,0.85,label) + #fig.text(0.1,0.85,label) + fig.text(0.535,0.12,label) + def annotateObs(): + fig = plt.gcf() + fig.text(0.1,0.85,r"Trenberth, K. E. and J. M. Caron, 2001: Estimates of Meridional Atmosphere and Ocean Heat Transports. J.Climate, 14, 3433-3443.", fontsize=8) + fig.text(0.1,0.825,r"Ganachaud, A. and C. Wunsch, 2000: Improved estimates of global ocean circulation, heat transport and mixing from hydrographic data.", fontsize=8) + fig.text(0.13,0.8,r"Nature, 408, 453-457", fontsize=8) + m6plot.setFigureSize(npanels=1) + # Load Observations + fObs = netCDF4.Dataset('/archive/John.Krasting/obs/TC2001/Trenberth_and_Caron_Heat_Transport.nc') #Trenberth and Caron + yobs = fObs.variables['ylat'][:] + NCEP = {}; NCEP['Global'] = fObs.variables['OTn'] + NCEP['Atlantic'] = fObs.variables['ATLn'][:]; NCEP['IndoPac'] = fObs.variables['INDPACn'][:] + ECMWF = {}; ECMWF['Global'] = fObs.variables['OTe'][:] + ECMWF['Atlantic'] = fObs.variables['ATLe'][:]; ECMWF['IndoPac'] = fObs.variables['INDPACe'][:] + + #G and W + Global = {} + Global['lat'] = numpy.array([-30., -19., 24., 47.]) + Global['trans'] = numpy.array([-0.6, -0.8, 1.8, 0.6]) + Global['err'] = numpy.array([0.3, 0.6, 0.3, 0.1]) + + Atlantic = {} + Atlantic['lat'] = numpy.array([-45., -30., -19., -11., -4.5, 7.5, 24., 47.]) + Atlantic['trans'] = numpy.array([0.66, 0.35, 0.77, 0.9, 1., 1.26, 1.27, 0.6]) + Atlantic['err'] = numpy.array([0.12, 0.15, 0.2, 0.4, 0.55, 0.31, 0.15, 0.09]) + + IndoPac = {} + IndoPac['lat'] = numpy.array([-30., -18., 24., 47.]) + IndoPac['trans'] = numpy.array([-0.9, -1.6, 0.52, 0.]) + IndoPac['err'] = numpy.array([0.3, 0.6, 0.2, 0.05,]) + + GandW = {} + GandW['Global'] = Global + GandW['Atlantic'] = Atlantic + GandW['IndoPac'] = IndoPac + + def plotGandW(lat,trans,err): + low = trans - err + high = trans + err + for n in range(0,len(low)): + if n == 0: + plt.plot([lat[n],lat[n]], [low[n],high[n]], 'c', linewidth=2.0, label='G&W') + else: + plt.plot([lat[n],lat[n]], [low[n],high[n]], 'c', linewidth=2.0) + plt.scatter(lat,trans,marker='s',facecolor='cyan') + # Global Heat Transport HTplot = heatTrans(advective,diffusive) yy = y[1:,:].max(axis=-1) plotHeatTrans(yy,HTplot,title='Global Y-Direction Heat Transport [PW]') + plt.plot(yobs,NCEP['Global'],'k--',linewidth=0.5,label='NCEP') + plt.plot(yobs,ECMWF['Global'],'k.',linewidth=0.5,label='ECMWF') + plotGandW(GandW['Global']['lat'],GandW['Global']['trans'],GandW['Global']['err']) plt.xlabel(r'Latitude [$\degree$N]') plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) + plt.legend(loc=0,fontsize=10) + annotateObs() if diffusive == None: annotatePlot('Warning: Diffusive component of transport is missing.') if stream != None: plt.savefig(stream[0]) @@ -106,8 +158,13 @@ def annotatePlot(label): yy = y[1:,:].max(axis=-1) HTplot[yy<-34] = numpy.nan plotHeatTrans(yy,HTplot,title='Atlantic Y-Direction Heat Transport [PW]') + plt.plot(yobs,NCEP['Atlantic'],'k--',linewidth=0.5,label='NCEP') + plt.plot(yobs,ECMWF['Atlantic'],'k.',linewidth=0.5,label='ECMWF') + plotGandW(GandW['Atlantic']['lat'],GandW['Atlantic']['trans'],GandW['Atlantic']['err']) plt.xlabel(r'Latitude [$\degree$N]') plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) + plt.legend(loc=0,fontsize=10) + annotateObs() if diffusive == None: annotatePlot('Warning: Diffusive component of transport is missing.') if stream != None: plt.savefig(stream[1]) @@ -119,11 +176,16 @@ def annotatePlot(label): m = 0*basin_code; m[(basin_code==3) | (basin_code==5)] = 1 HTplot = heatTrans(advective, diffusive, vmask=m*numpy.roll(m,-1,axis=-2)) yy = y[1:,:].max(axis=-1) - HTplot[yy<-38] = numpy.nan + HTplot[yy<-34] = numpy.nan plotHeatTrans(yy,HTplot,title='Indo-Pacific Y-Direction Heat Transport [PW]') + plt.plot(yobs,NCEP['IndoPac'],'k--',linewidth=0.5,label='NCEP') + plt.plot(yobs,ECMWF['IndoPac'],'k.',linewidth=0.5,label='ECMWF') + plotGandW(GandW['IndoPac']['lat'],GandW['IndoPac']['trans'],GandW['IndoPac']['err']) plt.xlabel(r'Latitude [$\degree$N]') + annotateObs() if diffusive == None: annotatePlot('Warning: Diffusive component of transport is missing.') plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) + plt.legend(loc=0,fontsize=10) if stream != None: plt.savefig(stream[2]) else: diff --git a/tools/analysis/zonal_S_bias_WOA05.py b/tools/analysis/zonal_S_bias_WOA05.py index a43c718576..0ac7f7f804 100755 --- a/tools/analysis/zonal_S_bias_WOA05.py +++ b/tools/analysis/zonal_S_bias_WOA05.py @@ -9,124 +9,143 @@ try: import argparse except: raise Exception('This version of python is not new enough. python 2.7 or newer is required.') -parser = argparse.ArgumentParser(description='''Script for plotting annual-average zonal salinity bias.''') -parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'salt' and 'e'.''') -parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') -parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') -parser.add_argument('-g','--gridspec', type=str, required=True, - help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') -parser.add_argument('-w','--woa', type=str, required=True, - help='''File containing WOA (or obs) data to compare against.''') -cmdLineArgs = parser.parse_args() - -if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') -if os.path.isdir(cmdLineArgs.gridspec): - xcenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] - y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2].max(axis=-1) - ycenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] - msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] - area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) - depth = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_topog.nc').variables['depth'][:] - try: basin = netCDF4.Dataset(cmdLineArgs.gridspec+'/basin_codes.nc').variables['basin'][:] - except: basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) -elif os.path.isfile(cmdLineArgs.gridspec): - xcenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[1::2,1::2] - y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2].max(axis=-1) - ycenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2] - msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] - area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) - depth = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_topog.nc','depth')[:] - try: basin = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'basin_codes.nc','basin')[:] - except: basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) -else: - raise ValueError('Unable to extract grid information from gridspec directory/tar file.') - -# Basin codes -basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) # All ocean points seeded from South Atlantic - -Sobs = netCDF4.Dataset( cmdLineArgs.woa ).variables['salt'] -if len(Sobs.shape)==3: Sobs = Sobs[:] -else: Sobs = Sobs[:].mean(axis=0) -Zobs = netCDF4.Dataset( cmdLineArgs.woa ).variables['eta'][:] - -rootGroup = netCDF4.Dataset( cmdLineArgs.annual_file ) -if 'salt' in rootGroup.variables: varName = 'salt' -elif 'so' in rootGroup.variables: varName = 'so' -else:raise Exception('Could not find "salt" or "so" in file "%s"'%(cmdLineArgs.annual_file)) -if len(rootGroup.variables[varName].shape)==4: Smod = rootGroup.variables[varName][:].mean(axis=0) -else: Smod = rootGroup.variables[varName][:] -if 'e' in rootGroup.variables: Zmod = rootGroup.variables['e'][0] -else: Zmod = Zobs # Using model z-output - -def zonalAverage(S, eta, area, mask=1.): - vols = ( mask * area ) * ( eta[:-1] - eta[1:] ) # mask * area * level thicknesses - return numpy.sum( vols * S, axis=-1 ) / numpy.sum( vols, axis=-1 ), (mask*eta).min(axis=-1) - -ci=m6plot.pmCI(0.125,2.25,.25) - -# Global -sPlot, z = zonalAverage(Smod, Zmod, area) -sObsPlot, _ = zonalAverage(Sobs, Zobs, area) -m6plot.yzplot( sPlot - sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title='''Global zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', - clim=ci, colormap='dunnePM', centerlabels=True, extend='both', - save=cmdLineArgs.outdir+'/S_global_xave_bias_WOA05.png') - -m6plot.yzcompare( sPlot, sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, - title1='Global zonal-average salinity [ppt]', - title2='''WOA'05 salinity [ppt]''', - clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', - dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, - save=cmdLineArgs.outdir+'/S_global_xave_bias_WOA05.3_panel.png') - -# Atlantic + Arctic -newMask = 1.*msk; newMask[ (basin!=2) & (basin!=4) ] = 0. -sPlot, z = zonalAverage(Smod, Zmod, area, mask=newMask) -sObsPlot, _ = zonalAverage(Sobs, Zobs, area, mask=newMask) -m6plot.yzplot( sPlot - sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title='''Atlantic zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', - clim=ci, colormap='dunnePM', centerlabels=True, extend='both', - save=cmdLineArgs.outdir+'/S_Atlantic_xave_bias_WOA05.png') - -m6plot.yzcompare( sPlot, sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, - title1='Atlantic zonal-average salinity [ppt]', - title2='''WOA'05 salinity [ppt]''', - clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', - dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, - save=cmdLineArgs.outdir+'/S_Atlantic_xave_bias_WOA05.3_panel.png') - -# Pacific -newMask = 1.*msk; newMask[ (basin!=3) ] = 0. -sPlot, z = zonalAverage(Smod, Zmod, area, mask=newMask) -sObsPlot, _ = zonalAverage(Sobs, Zobs, area, mask=newMask) -m6plot.yzplot( sPlot - sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title='''Pacific zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', - clim=ci, colormap='dunnePM', centerlabels=True, extend='both', - save=cmdLineArgs.outdir+'/S_Pacific_xave_bias_WOA05.png') - -m6plot.yzcompare( sPlot, sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, - title1='Pacific zonal-average salinity [ppt]', - title2='''WOA'05 salinity [ppt]''', - clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', - dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, - save=cmdLineArgs.outdir+'/S_Pacific_xave_bias_WOA05.3_panel.png') - -# Indian -newMask = 1.*msk; newMask[ (basin!=5) ] = 0. -sPlot, z = zonalAverage(Smod, Zmod, area, mask=newMask) -sObsPlot, _ = zonalAverage(Sobs, Zobs, area, mask=newMask) -m6plot.yzplot( sPlot - sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title='''Indian zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', - clim=ci, colormap='dunnePM', centerlabels=True, extend='both', - save=cmdLineArgs.outdir+'/S_Indian_xave_bias_WOA05.png') - -m6plot.yzcompare( sPlot, sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, - title1='Indian zonal-average salinity [ppt]', - title2='''WOA'05 salinity [ppt]''', - clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', - dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, - save=cmdLineArgs.outdir+'/S_Indian_xave_bias_WOA05.3_panel.png') +def run(): + parser = argparse.ArgumentParser(description='''Script for plotting annual-average zonal salinity bias.''') + parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'salt' and 'e'.''') + parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') + parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') + parser.add_argument('-g','--gridspec', type=str, required=True, + help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') + parser.add_argument('-w','--woa', type=str, required=True, + help='''File containing WOA (or obs) data to compare against.''') + cmdLineArgs = parser.parse_args() + main(cmdLineArgs) + +def main(cmdLineArgs,stream=None): + if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') + if os.path.isdir(cmdLineArgs.gridspec): + xcenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] + y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2].max(axis=-1) + ycenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] + msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] + area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_topog.nc').variables['depth'][:] + try: basin = netCDF4.Dataset(cmdLineArgs.gridspec+'/basin_codes.nc').variables['basin'][:] + except: basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) + elif os.path.isfile(cmdLineArgs.gridspec): + xcenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[1::2,1::2] + y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2].max(axis=-1) + ycenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2] + msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] + area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_topog.nc','depth')[:] + try: basin = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'basin_codes.nc','basin')[:] + except: basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) + else: + raise ValueError('Unable to extract grid information from gridspec directory/tar file.') + + if stream != None: + if len(stream) != 4: + raise ValueError('If specifying output streams, exactly four (4) streams are needed for this analysis') + + Sobs = netCDF4.Dataset( cmdLineArgs.woa ).variables['salt'] + if len(Sobs.shape)==3: Sobs = Sobs[:] + else: Sobs = Sobs[:].mean(axis=0) + Zobs = netCDF4.Dataset( cmdLineArgs.woa ).variables['eta'][:] + + rootGroup = netCDF4.MFDataset( cmdLineArgs.annual_file ) + if 'salt' in rootGroup.variables: varName = 'salt' + elif 'so' in rootGroup.variables: varName = 'so' + else:raise Exception('Could not find "salt" or "so" in file "%s"'%(cmdLineArgs.annual_file)) + if len(rootGroup.variables[varName].shape)==4: Smod = rootGroup.variables[varName][:].mean(axis=0) + else: Smod = rootGroup.variables[varName][:] + if 'e' in rootGroup.variables: Zmod = rootGroup.variables['e'][0] + else: Zmod = Zobs # Using model z-output + + def zonalAverage(S, eta, area, mask=1.): + vols = ( mask * area ) * ( eta[:-1] - eta[1:] ) # mask * area * level thicknesses + return numpy.sum( vols * S, axis=-1 ) / numpy.sum( vols, axis=-1 ), (mask*eta).min(axis=-1) + + ci=m6plot.pmCI(0.125,2.25,.25) + + # Global + sPlot, z = zonalAverage(Smod, Zmod, area) + sObsPlot, _ = zonalAverage(Sobs, Zobs, area) + if stream != None: objOut = stream[0] + else: objOut = cmdLineArgs.outdir+'/S_global_xave_bias_WOA05.png' + m6plot.yzplot( sPlot - sObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, title='''Global zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', + clim=ci, colormap='dunnePM', centerlabels=True, extend='both', + save=objOut) + + if stream == None: + m6plot.yzcompare( sPlot, sObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, + title1='Global zonal-average salinity [ppt]', + title2='''WOA'05 salinity [ppt]''', + clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', + dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, + save=cmdLineArgs.outdir+'/S_global_xave_bias_WOA05.3_panel.png') + + # Atlantic + Arctic + newMask = 1.*msk; newMask[ (basin!=2) & (basin!=4) ] = 0. + sPlot, z = zonalAverage(Smod, Zmod, area, mask=newMask) + sObsPlot, _ = zonalAverage(Sobs, Zobs, area, mask=newMask) + if stream != None: objOut = stream[1] + else: objOut = cmdLineArgs.outdir+'/S_Atlantic_xave_bias_WOA05.png' + m6plot.yzplot( sPlot - sObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, title='''Atlantic zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', + clim=ci, colormap='dunnePM', centerlabels=True, extend='both', + save=objOut) + + if stream == None: + m6plot.yzcompare( sPlot, sObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, + title1='Atlantic zonal-average salinity [ppt]', + title2='''WOA'05 salinity [ppt]''', + clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', + dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, + save=cmdLineArgs.outdir+'/S_Atlantic_xave_bias_WOA05.3_panel.png') + + # Pacific + newMask = 1.*msk; newMask[ (basin!=3) ] = 0. + sPlot, z = zonalAverage(Smod, Zmod, area, mask=newMask) + sObsPlot, _ = zonalAverage(Sobs, Zobs, area, mask=newMask) + if stream != None: objOut = stream[2] + else: objOut = cmdLineArgs.outdir+'/S_Pacific_xave_bias_WOA05.png' + m6plot.yzplot( sPlot - sObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, title='''Pacific zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', + clim=ci, colormap='dunnePM', centerlabels=True, extend='both', + save=objOut) + + if stream == None: + m6plot.yzcompare( sPlot, sObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, + title1='Pacific zonal-average salinity [ppt]', + title2='''WOA'05 salinity [ppt]''', + clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', + dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, + save=cmdLineArgs.outdir+'/S_Pacific_xave_bias_WOA05.3_panel.png') + + # Indian + newMask = 1.*msk; newMask[ (basin!=5) ] = 0. + sPlot, z = zonalAverage(Smod, Zmod, area, mask=newMask) + sObsPlot, _ = zonalAverage(Sobs, Zobs, area, mask=newMask) + if stream != None: objOut = stream[3] + else: objOut = cmdLineArgs.outdir+'/S_Indian_xave_bias_WOA05.png' + m6plot.yzplot( sPlot - sObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, title='''Indian zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', + clim=ci, colormap='dunnePM', centerlabels=True, extend='both', + save=objOut) + + if stream == None: + m6plot.yzcompare( sPlot, sObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, + title1='Indian zonal-average salinity [ppt]', + title2='''WOA'05 salinity [ppt]''', + clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', + dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, + save=cmdLineArgs.outdir+'/S_Indian_xave_bias_WOA05.3_panel.png') + +if __name__ == '__main__': + run() diff --git a/tools/analysis/zonal_T_bias_WOA05.py b/tools/analysis/zonal_T_bias_WOA05.py index 6b7b36fc3f..66c36540ea 100755 --- a/tools/analysis/zonal_T_bias_WOA05.py +++ b/tools/analysis/zonal_T_bias_WOA05.py @@ -9,124 +9,146 @@ try: import argparse except: raise Exception('This version of python is not new enough. python 2.7 or newer is required.') -parser = argparse.ArgumentParser(description='''Script for plotting annual-average zonal temperature bias.''') -parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'temp' and 'e'.''') -parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') -parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') -parser.add_argument('-g','--gridspec', type=str, required=True, - help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') -parser.add_argument('-w','--woa', type=str, required=True, - help='''File containing WOA (or obs) data to compare against.''') -cmdLineArgs = parser.parse_args() - -if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') -if os.path.isdir(cmdLineArgs.gridspec): - xcenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] - y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2].max(axis=-1) - ycenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] - msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] - area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) - depth = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_topog.nc').variables['depth'][:] - try: basin = netCDF4.Dataset(cmdLineArgs.gridspec+'/basin_codes.nc').variables['basin'][:] - except: basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) -elif os.path.isfile(cmdLineArgs.gridspec): - xcenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[1::2,1::2] - y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2].max(axis=-1) - ycenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2] - msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] - area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) - depth = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_topog.nc','depth')[:] - try: basin = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'basin_codes.nc','basin')[:] - except: basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) -else: - raise ValueError('Unable to extract grid information from gridspec directory/tar file.') - -Tobs = netCDF4.Dataset( cmdLineArgs.woa ) -if 'temp' in Tobs.variables: Tobs = Tobs.variables['temp'] -else: Tobs = Tobs.variables['ptemp'] -if len(Tobs.shape)==3: Tobs = Tobs[:] -else: Tobs = Tobs[:].mean(axis=0) -Zobs = netCDF4.Dataset( cmdLineArgs.woa ).variables['eta'][:] - -rootGroup = netCDF4.Dataset( cmdLineArgs.annual_file ) -if 'temp' in rootGroup.variables: varName = 'temp' -elif 'ptemp' in rootGroup.variables: varName = 'ptemp' -elif 'thetao' in rootGroup.variables: varName = 'thetao' -else: raise Exception('Could not find "temp", "ptemp" or "thetao" in file "%s"'%(cmdLineArgs.annual_file)) -if len(rootGroup.variables[varName].shape)==4: Tmod = rootGroup.variables[varName][:].mean(axis=0) -else: Tmod = rootGroup.variables[varName][:] -if 'e' in rootGroup.variables: Zmod = rootGroup.variables['e'][0] -else: Zmod = Zobs # Using model z-ou:put - -def zonalAverage(T, eta, area, mask=1.): - vols = ( mask * area ) * ( eta[:-1] - eta[1:] ) # mask * area * level thicknesses - return numpy.sum( vols * T, axis=-1 ) / numpy.sum( vols, axis=-1 ), (mask*eta).min(axis=-1) - -ci=m6plot.pmCI(0.25,4.5,.5) - -# Global -tPlot, z = zonalAverage(Tmod, Zmod, area) -tObsPlot, _ = zonalAverage(Tobs, Zobs, area) -m6plot.yzplot( tPlot - tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title=r'''Global zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', - clim=ci, colormap='dunnePM', centerlabels=True, extend='both', - save=cmdLineArgs.outdir+'/T_global_xave_bias_WOA05.png') - -m6plot.yzcompare( tPlot, tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, - title1=r'Global zonal-average $\theta$ [$\degree$C]', - title2=r'''WOA'05 $\theta$ [$\degree$C]''', - clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', - dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, - save=cmdLineArgs.outdir+'/T_global_xave_bias_WOA05.3_panel.png') - -# Atlantic + Arctic -newMask = 1.*msk; newMask[ (basin!=2) & (basin!=4) ] = 0. -tPlot, z = zonalAverage(Tmod, Zmod, area, mask=newMask) -tObsPlot, _ = zonalAverage(Tobs, Zobs, area, mask=newMask) -m6plot.yzplot( tPlot - tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title=r'''Atlantic zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', - clim=ci, colormap='dunnePM', centerlabels=True, extend='both', - save=cmdLineArgs.outdir+'/T_Atlantic_xave_bias_WOA05.png') - -m6plot.yzcompare( tPlot, tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, - title1=r'Atlantic zonal-average $\theta$ [$\degree$C]', - title2=r'''WOA'05 $\theta$ [$\degree$C]''', - clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', - dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, - save=cmdLineArgs.outdir+'/T_Atlantic_xave_bias_WOA05.3_panel.png') - -# Pacific -newMask = 1.*msk; newMask[ (basin!=3) ] = 0. -tPlot, z = zonalAverage(Tmod, Zmod, area, mask=newMask) -tObsPlot, _ = zonalAverage(Tobs, Zobs, area, mask=newMask) -m6plot.yzplot( tPlot - tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title=r'''Pacific zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', - clim=ci, colormap='dunnePM', centerlabels=True, extend='both', - save=cmdLineArgs.outdir+'/T_Pacific_xave_bias_WOA05.png') - -m6plot.yzcompare( tPlot, tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, - title1=r'Pacific zonal-average $\theta$ [$\degree$C]', - title2=r'''WOA'05 $\theta$ [$\degree$C]''', - clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', - dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, - save=cmdLineArgs.outdir+'/T_Pacific_xave_bias_WOA05.3_panel.png') - -# Indian -newMask = 1.*msk; newMask[ (basin!=5) ] = 0. -tPlot, z = zonalAverage(Tmod, Zmod, area, mask=newMask) -tObsPlot, _ = zonalAverage(Tobs, Zobs, area, mask=newMask) -m6plot.yzplot( tPlot - tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title=r'''Indian zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', - clim=ci, colormap='dunnePM', centerlabels=True, extend='both', - save=cmdLineArgs.outdir+'/T_Indian_xave_bias_WOA05.png') - -m6plot.yzcompare( tPlot, tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, - title1=r'Indian zonal-average $\theta$ [$\degree$C]', - title2=r'''WOA'05 $\theta$ [$\degree$C]''', - clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', - dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, - save=cmdLineArgs.outdir+'/T_Indian_xave_bias_WOA05.3_panel.png') +def run(): + parser = argparse.ArgumentParser(description='''Script for plotting annual-average zonal temperature bias.''') + parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'temp' and 'e'.''') + parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') + parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') + parser.add_argument('-g','--gridspec', type=str, required=True, + help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') + parser.add_argument('-w','--woa', type=str, required=True, + help='''File containing WOA (or obs) data to compare against.''') + cmdLineArgs = parser.parse_args() + main(cmdLineArgs) + +def main(cmdLineArgs,stream=None): + if not os.path.exists(cmdLineArgs.gridspec): raise ValueError('Specified gridspec directory/tar file does not exist.') + if os.path.isdir(cmdLineArgs.gridspec): + xcenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][1::2,1::2] + y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2].max(axis=-1) + ycenter = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][1::2,1::2] + msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] + area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_topog.nc').variables['depth'][:] + try: basin = netCDF4.Dataset(cmdLineArgs.gridspec+'/basin_codes.nc').variables['basin'][:] + except: basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) + elif os.path.isfile(cmdLineArgs.gridspec): + xcenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','x')[1::2,1::2] + y = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2].max(axis=-1) + ycenter = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','y')[1::2,1::2] + msk = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_mask.nc','mask')[:] + area = msk*m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_hgrid.nc','area')[:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) + depth = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'ocean_topog.nc','depth')[:] + try: basin = m6toolbox.readNCFromTar(cmdLineArgs.gridspec,'basin_codes.nc','basin')[:] + except: basin = m6toolbox.genBasinMasks(xcenter, ycenter, depth) + else: + raise ValueError('Unable to extract grid information from gridspec directory/tar file.') + + if stream != None: + if len(stream) != 4: + raise ValueError('If specifying output streams, exactly four (4) streams are needed for this analysis') + + Tobs = netCDF4.Dataset( cmdLineArgs.woa ) + if 'temp' in Tobs.variables: Tobs = Tobs.variables['temp'] + else: Tobs = Tobs.variables['ptemp'] + if len(Tobs.shape)==3: Tobs = Tobs[:] + else: Tobs = Tobs[:].mean(axis=0) + Zobs = netCDF4.Dataset( cmdLineArgs.woa ).variables['eta'][:] + + rootGroup = netCDF4.MFDataset( cmdLineArgs.annual_file ) + if 'temp' in rootGroup.variables: varName = 'temp' + elif 'ptemp' in rootGroup.variables: varName = 'ptemp' + elif 'thetao' in rootGroup.variables: varName = 'thetao' + else: raise Exception('Could not find "temp", "ptemp" or "thetao" in file "%s"'%(cmdLineArgs.annual_file)) + if len(rootGroup.variables[varName].shape)==4: Tmod = rootGroup.variables[varName][:].mean(axis=0) + else: Tmod = rootGroup.variables[varName][:] + if 'e' in rootGroup.variables: Zmod = rootGroup.variables['e'][0] + else: Zmod = Zobs # Using model z-ou:put + + def zonalAverage(T, eta, area, mask=1.): + vols = ( mask * area ) * ( eta[:-1] - eta[1:] ) # mask * area * level thicknesses + return numpy.sum( vols * T, axis=-1 ) / numpy.sum( vols, axis=-1 ), (mask*eta).min(axis=-1) + + ci=m6plot.pmCI(0.25,4.5,.5) + + # Global + tPlot, z = zonalAverage(Tmod, Zmod, area) + tObsPlot, _ = zonalAverage(Tobs, Zobs, area) + if stream != None: objOut = stream[0] + else: objOut = cmdLineArgs.outdir+'/T_global_xave_bias_WOA05.png' + m6plot.yzplot( tPlot - tObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, title=r'''Global zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', + clim=ci, colormap='dunnePM', centerlabels=True, extend='both', + save=objOut) + + if stream == None: + m6plot.yzcompare( tPlot, tObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, + title1=r'Global zonal-average $\theta$ [$\degree$C]', + title2=r'''WOA'05 $\theta$ [$\degree$C]''', + clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', + dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, + save=cmdLineArgs.outdir+'/T_global_xave_bias_WOA05.3_panel.png') + + # Atlantic + Arctic + newMask = 1.*msk; newMask[ (basin!=2) & (basin!=4) ] = 0. + tPlot, z = zonalAverage(Tmod, Zmod, area, mask=newMask) + tObsPlot, _ = zonalAverage(Tobs, Zobs, area, mask=newMask) + if stream != None: objOut = stream[1] + else: objOut = cmdLineArgs.outdir+'/T_Atlantic_xave_bias_WOA05.png' + m6plot.yzplot( tPlot - tObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, title=r'''Atlantic zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', + clim=ci, colormap='dunnePM', centerlabels=True, extend='both', + save=objOut) + + if stream == None: + m6plot.yzcompare( tPlot, tObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, + title1=r'Atlantic zonal-average $\theta$ [$\degree$C]', + title2=r'''WOA'05 $\theta$ [$\degree$C]''', + clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', + dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, + save=cmdLineArgs.outdir+'/T_Atlantic_xave_bias_WOA05.3_panel.png') + + # Pacific + newMask = 1.*msk; newMask[ (basin!=3) ] = 0. + tPlot, z = zonalAverage(Tmod, Zmod, area, mask=newMask) + tObsPlot, _ = zonalAverage(Tobs, Zobs, area, mask=newMask) + if stream != None: objOut = stream[2] + else: objOut = cmdLineArgs.outdir+'/T_Pacific_xave_bias_WOA05.png' + m6plot.yzplot( tPlot - tObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, title=r'''Pacific zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', + clim=ci, colormap='dunnePM', centerlabels=True, extend='both', + save=objOut) + + if stream == None: + m6plot.yzcompare( tPlot, tObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, + title1=r'Pacific zonal-average $\theta$ [$\degree$C]', + title2=r'''WOA'05 $\theta$ [$\degree$C]''', + clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', + dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, + save=cmdLineArgs.outdir+'/T_Pacific_xave_bias_WOA05.3_panel.png') + + # Indian + newMask = 1.*msk; newMask[ (basin!=5) ] = 0. + tPlot, z = zonalAverage(Tmod, Zmod, area, mask=newMask) + tObsPlot, _ = zonalAverage(Tobs, Zobs, area, mask=newMask) + if stream != None: objOut = stream[3] + else: objOut = cmdLineArgs.outdir+'/T_Indian_xave_bias_WOA05.png' + m6plot.yzplot( tPlot - tObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, title=r'''Indian zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', + clim=ci, colormap='dunnePM', centerlabels=True, extend='both', + save=objOut) + + if stream == None: + m6plot.yzcompare( tPlot, tObsPlot , y, z, splitscale=[0., -1000., -6500.], + suptitle=rootGroup.title+' '+cmdLineArgs.label, + title1=r'Indian zonal-average $\theta$ [$\degree$C]', + title2=r'''WOA'05 $\theta$ [$\degree$C]''', + clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', + dlim=ci, dcolormap='dunnePM', dextend='both', centerdlabels=True, + save=cmdLineArgs.outdir+'/T_Indian_xave_bias_WOA05.3_panel.png') + +if __name__ == '__main__': + run() From bd77926dc4a72cb8d6fad9e3c22a72c0925a59cb Mon Sep 17 00:00:00 2001 From: John Krasting Date: Tue, 1 Mar 2016 09:58:51 -0500 Subject: [PATCH 7/9] Simple CSH script to repair automated frepp analysis - Intended to facilitate regeneration of analysis figures for earlier 0.5 degree simulations --- tools/analysis/manual_frepp_analysis.csh | 69 ++++++++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100755 tools/analysis/manual_frepp_analysis.csh diff --git a/tools/analysis/manual_frepp_analysis.csh b/tools/analysis/manual_frepp_analysis.csh new file mode 100755 index 0000000000..08b0a57fc4 --- /dev/null +++ b/tools/analysis/manual_frepp_analysis.csh @@ -0,0 +1,69 @@ +#!/bin/csh -f + +# ------------------------------------------------- +# Simple csh snippet to repair frepp analysis +# +# *** make sure to run 'module load python' *** +# +# Instructions: +# 1. Fill in ppDir, years, and chunk +# 2. Select either the 0.25 or 0.5 degree +# GRIDSPEC and WOA paths +# +# ------------------------------------------------- + + +set ppDir = "" # e.g. include everything up to but not including the '/pp' +set yr1 = "" # make sure it has 4 digits +set yr2 = "" # make sure it has 4 digits +set chunk = "" # e.g. '${chunk}' or '20yr' + +# Settings for the 0.5 degree configuration +set GRIDSPEC = "/archive/gold/datasets/OM4_05/mosaic.v20151203.unpacked" +set WOA = "/archive/gold/datasets/OM4_05/obs/WOA05_ptemp_salt_annual.v2015.12.03.nc" + +# Settings for the 0.25 degree configuration +# set GRIDSPEC = "/archive/gold/datasets/OM4_025/mosaic.v20140610.tar" +# set WOA = "/archive/gold/datasets/OM4_025/obs/WOA05_ptemp_salt_annual.v20150310.nc" + + +# No need to edit below this line +# ----------------------------------------------------------------------------- + +mkdir -p ../../../ocean_${yr1}-${yr2}/ptemp +./SST_bias_WOA05.py -w ${WOA} -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual_z/av/annual_${chunk}/ocean_annual_z.${yr1}-${yr2}.ann.nc +./zonal_T_bias_WOA05.py -w ${WOA} -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual_z/av/annual_${chunk}/ocean_annual_z.${yr1}-${yr2}.ann.nc +./depth_average_T_bias.py -zb 100 -w ${WOA} -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual_z/av/annual_${chunk}/ocean_annual_z.${yr1}-${yr2}.ann.nc +./depth_average_T_bias.py -zb 300 -w ${WOA} -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual_z/av/annual_${chunk}/ocean_annual_z.${yr1}-${yr2}.ann.nc +./depth_average_T_bias.py -zb 700 -w ${WOA} -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual_z/av/annual_${chunk}/ocean_annual_z.${yr1}-${yr2}.ann.nc +./depth_average_T_bias.py -zb 2000 -w ${WOA} -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual_z/av/annual_${chunk}/ocean_annual_z.${yr1}-${yr2}.ann.nc +./depth_average_T_bias.py -zt 2000 -zb 4000 -w ${WOA} -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual_z/av/annual_${chunk}/ocean_annual_z.${yr1}-${yr2}.ann.nc +./depth_average_T_bias.py -zb 6500 -w ${WOA} -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/ptemp -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual_z/av/annual_${chunk}/ocean_annual_z.${yr1}-${yr2}.ann.nc + + +mkdir -p ../../../ocean_${yr1}-${yr2}/salinity +./SSS_bias_WOA05.py -w ${WOA} -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/salinity -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual_z/av/annual_${chunk}/ocean_annual_z.${yr1}-${yr2}.ann.nc +./zonal_S_bias_WOA05.py -w ${WOA} -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/salinity -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual_z/av/annual_${chunk}/ocean_annual_z.${yr1}-${yr2}.ann.nc + + +# -- This script may not work; different behavior on workstation vs. PAN +mkdir -p ../../../ocean_${yr1}-${yr2}/sections +./vertical_sections_annual_bias_WOA05.py -w ${WOA} -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/sections -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual_z/av/annual_${chunk}/ocean_annual_z.${yr1}-${yr2}.ann.nc + + +mkdir -p ../../../ocean_${yr1}-${yr2}/MOC +./meridional_overturning.py -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/MOC -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual_z/av/annual_${chunk}/ocean_annual_z.${yr1}-${yr2}.ann.nc + + +mkdir -p ../../../ocean_${yr1}-${yr2}/Hosoda_MLD +./MLD_003.py -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/Hosoda_MLD -l ${yr1}-${yr2} ${ppDir}/pp/ocean_monthly/ts/monthly/${chunk}/ocean_monthly.${yr1}01-${yr2}12.MLD_003.nc + + +mkdir -p ../../../ocean_${yr1}-${yr2}/heat_transport +./poleward_heat_transport.py -l ${yr1}-${yr2} -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/heat_transport -l ${yr1}-${yr2} ${ppDir}/pp/ocean_annual/av/annual_${chunk}/ocean_annual.${yr1}-${yr2}.ann.nc + +# -- This script may not work; requires MIDAS build +mkdir -p ../../../ocean_${yr1}-${yr2}/heat_salt_0_300m +./TS_depth_integrals.py -r ${WOA} -s 0 -e 300 -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/heat_salt_0_300m -l ${yr1}-${yr2} ${ppDir}/ocean_annual.$yr1-$yr2.ann.nc + +#// From 0b71b84ee4430dfaec175e5f4f2209462a1941b9 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Wed, 2 Mar 2016 08:13:54 -0500 Subject: [PATCH 8/9] Set-up modules environment for GFDL systems --- tools/analysis/manual_frepp_analysis.csh | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tools/analysis/manual_frepp_analysis.csh b/tools/analysis/manual_frepp_analysis.csh index 08b0a57fc4..4528737798 100755 --- a/tools/analysis/manual_frepp_analysis.csh +++ b/tools/analysis/manual_frepp_analysis.csh @@ -1,9 +1,8 @@ -#!/bin/csh -f +#!/bin/csh # ------------------------------------------------- # Simple csh snippet to repair frepp analysis # -# *** make sure to run 'module load python' *** # # Instructions: # 1. Fill in ppDir, years, and chunk @@ -12,6 +11,8 @@ # # ------------------------------------------------- +source $MODULESHOME/init/csh +module load python set ppDir = "" # e.g. include everything up to but not including the '/pp' set yr1 = "" # make sure it has 4 digits From 08499ae11f15d730e5188ef74426d97aa73af418 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Wed, 2 Mar 2016 08:22:30 -0500 Subject: [PATCH 9/9] Added error messages for some common failures --- tools/analysis/manual_frepp_analysis.csh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tools/analysis/manual_frepp_analysis.csh b/tools/analysis/manual_frepp_analysis.csh index 4528737798..36707e6da3 100755 --- a/tools/analysis/manual_frepp_analysis.csh +++ b/tools/analysis/manual_frepp_analysis.csh @@ -58,6 +58,7 @@ mkdir -p ../../../ocean_${yr1}-${yr2}/MOC mkdir -p ../../../ocean_${yr1}-${yr2}/Hosoda_MLD ./MLD_003.py -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/Hosoda_MLD -l ${yr1}-${yr2} ${ppDir}/pp/ocean_monthly/ts/monthly/${chunk}/ocean_monthly.${yr1}01-${yr2}12.MLD_003.nc +if ($status != 0) echo ' --> Mixed Layer Depth analysis may only be available on 5-year chunks.' mkdir -p ../../../ocean_${yr1}-${yr2}/heat_transport @@ -66,5 +67,6 @@ mkdir -p ../../../ocean_${yr1}-${yr2}/heat_transport # -- This script may not work; requires MIDAS build mkdir -p ../../../ocean_${yr1}-${yr2}/heat_salt_0_300m ./TS_depth_integrals.py -r ${WOA} -s 0 -e 300 -g ${GRIDSPEC} -o ../../../ocean_${yr1}-${yr2}/heat_salt_0_300m -l ${yr1}-${yr2} ${ppDir}/ocean_annual.$yr1-$yr2.ann.nc +if ($status != 0) echo ' --> Script may have failed on the MIDAS dependency.' #//