From 47cb72e7a644c415ccc3e121a01c0b29c2d659d4 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Fri, 15 Sep 2017 15:33:45 -0400 Subject: [PATCH 01/22] Initial commit of CMIP6 offline diag script - basin average overturning - basin average heat transports --- tools/analysis/CMIP6_offline_diags.py | 235 ++++++++++++++++++++++++++ tools/analysis/m6toolbox.py | 28 +++ 2 files changed, 263 insertions(+) create mode 100755 tools/analysis/CMIP6_offline_diags.py diff --git a/tools/analysis/CMIP6_offline_diags.py b/tools/analysis/CMIP6_offline_diags.py new file mode 100755 index 0000000000..5039c00cce --- /dev/null +++ b/tools/analysis/CMIP6_offline_diags.py @@ -0,0 +1,235 @@ +#!/usr/bin/env python + +import argparse +import m6toolbox +import netCDF4 as nc +import numpy as np +import os +import sys +import matplotlib.pyplot as plt + +##-- RefineDiag Script for CMIP6 +## +## Variables we intend to provide in z-coordinates: +## +## msftyyz -> vmo (ocean_z) * both 0.25 and 0.5 resolutions +## msftyzmpa -> vhGM (ocean_z) * applies only to 0.5 resolution +## msftyzsmpa -> vhml (ocean_z) * both 0.25 and 0.5 resolutions +## +## +## Variables we intend to provide in rho-coordinates: +## (potenital density referenced to 2000 m) +## +## msftyrho -> vmo +## msftyrhompa -> vhGM * applies only to 0.5 resolution +## +## +## 2-D variables we intent to provide: +## +## hfy -> T_ady_2d + ndiff_tracer_trans_y_2d_T * T_ady_2d needs to be converted to Watts (Cp = 3992.) +## ndiff_tracer_trans_y_2d_T already in Watts +## advective term in both 0.25 and 0.5 resolutions +## diffusive term only in 0.5 resolution +## +## hfx -> same recipie as above, expect for x-dimension +## hfbasin -> summed line of hfy in each basin +## +## +## Outstanding issues +## 1.) regirdding of vh, vhGM to rho-corrdinates +## 2.) vhGM and vhML units need to be in kg s-1 +## 2.) save out RHO_0 and Cp somewhere in netCDF files to key off of +## +## +## CMIP variables that will NOT be provided: +## +## hfbasinpmadv, hfbasinpsmadv, hfbasinpmdiff, hfbasinpadv +## (We advect the tracer with the residual mean velocity; individual terms cannot be diagnosed) +## +## htovgyre, htovovrt, sltovgyre, sltovovrt +## (Code for offline calculation not ready.) +## +##-- + +def run(): + parser = argparse.ArgumentParser(description='''CMIP6 RefineDiag Script for OM4''') + parser.add_argument('infile', type=str, help='''Input file''') + parser.add_argument('-b','--basinfile', type=str, default='', required=True, help='''File containing OM4 basin masks''') + parser.add_argument('-o','--outfile', type=str, default=None, help='''Output file name''') + args = parser.parse_args() + main(args) + +def main(args): + #-- Define Regions and their associated masks + # Note: The Atlantic should include other smaller bays/seas that are + # included in the definition used in meridional_overturning.py + + region = np.array(['atlantic_arctic_ocean','indian_pacific_ocean','global_ocean']) + + #-- Read basin masks + f_basin = nc.Dataset(args.basinfile) + basin_code = f_basin.variables['basin'][:] + + atlantic_arctic_mask = basin_code * 0. + atlantic_arctic_mask[(basin_code==2) | (basin_code==4) | (basin_code==6) | (basin_code==7) | (basin_code==8)] = 1. + + indo_pacific_mask = basin_code * 0. + indo_pacific_mask[(basin_code==3) | (basin_code==5)] = 1. + + #-- Read model data + f_in = nc.Dataset(args.infile) + + #-- Read in existing dimensions from history netcdf file + yh = f_in.variables['yh'] + z_l = f_in.variables['z_l'] + z_i = f_in.variables['z_i'] + tax = f_in.variables['time'] + + #-- msftyyz + varname = 'vmo' + msftyyz = np.ma.ones((len(tax),3,len(z_l),len(yh)))*0. + msftyyz[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) + msftyyz[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) + msftyyz[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) + msftyyz[msftyyz.mask] = 1.e20 + msftyyz = np.ma.array(msftyyz,fill_value=1.e20) + msftyyz.long_name = 'Ocean Y Overturning Mass Streamfunction' + msftyyz.units = 'kg s-1' + msftyyz.cell_methods = 'z_l:sum yh:sum basin:mean time:mean' + msftyyz.time_avg_info = 'average_T1,average_T2,average_DT' + msftyyz.standard_name = 'ocean_y_overturning_mass_streamfunction' + + #-- msftyzmpa + varname = 'vhGM' + msftyzmpa = np.ma.ones((len(tax),3,len(z_l),len(yh)))*0. + msftyzmpa[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) + msftyzmpa[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) + msftyzmpa[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) + msftyzmpa[msftyzmpa.mask] = 1.e20 + msftyzmpa = np.ma.array(msftyzmpa,fill_value=1.e20) + msftyzmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized mesoscale advection' + msftyzmpa.units = 'kg s-1' + msftyzmpa.cell_methods = 'z_l:sum yh:sum basin:mean time:mean' + msftyzmpa.time_avg_info = 'average_T1,average_T2,average_DT' + msftyzmpa.standard_name = 'ocean_meridional_overturning_mass_streamfunction_due_to_parameterized_'+\ + 'mesoscale_advection' + + #-- msftyzsmpa + varname = 'vhml' + msftyzsmpa = np.ma.ones((len(tax),3,len(z_l),len(yh)))*0. + msftyzsmpa[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) + msftyzsmpa[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) + msftyzsmpa[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) + msftyzsmpa[msftyzsmpa.mask] = 1.e20 + msftyzsmpa = np.ma.array(msftyzsmpa,fill_value=1.e20) + msftyzsmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized submesoscale advection' + msftyzsmpa.units = 'kg s-1' + msftyzsmpa.cell_methods = 'z_l:sum yh:sum basin:mean time:mean' + msftyzsmpa.time_avg_info = 'average_T1,average_T2,average_DT' + msftyzsmpa.standard_name = 'ocean_meridional_overturning_mass_streamfunction_due_to_parameterized_'+\ + 'submesoscale_advection' + + #-- Read time bounds + nv = f_in.variables['nv'] + average_T1 = f_in.variables['average_T1'] + average_T2 = f_in.variables['average_T2'] + average_DT = f_in.variables['average_DT'] + time_bnds = f_in.variables['time_bnds'] + + #-- Generate output filename + if args.outfile is None: + if hasattr(f_in,'filename'): + args.outfile = f_in.filename + else: + args.outfile = os.path.basename(args.infile) + args.outfile = args.outfile.split('.') + args.outfile[-2] = args.outfile[-2]+'_refined' + args.outfile = '.'.join(args.outfile) + + #-- Write output file + try: + os.remove(args.outfile) + except: + pass + + if os.path.exists(args.outfile): + raise IOError('Output netCDF file already exists.') + exit(1) + + f_out = nc.Dataset(args.outfile, 'w', format='NETCDF3_CLASSIC') + f_out.setncatts(f_in.__dict__) + f_out.filename = args.outfile + + time_dim = f_out.createDimension('time', size=None) + basin_dim = f_out.createDimension('basin', size=3) + strlen_dim = f_out.createDimension('strlen', size=21) + yh_dim = f_out.createDimension('yh', size=len(yh[:])) + z_l_dim = f_out.createDimension('z_l', size=len(z_l[:])) + z_i_dim = f_out.createDimension('z_i', size=len(z_i[:])) + nv_dim = f_out.createDimension('nv', size=len(nv[:])) + + time_out = f_out.createVariable('time', np.float64, ('time')) + #basin_out = f_out.createVariable('basin', np.int32, ('basin')) + yh_out = f_out.createVariable('yh', np.float64, ('yh')) + region_out = f_out.createVariable('region', 'c', ('basin', 'strlen')) + z_l_out = f_out.createVariable('z_l', np.float64, ('z_l')) + z_i_out = f_out.createVariable('z_i', np.float64, ('z_i')) + nv_out = f_out.createVariable('nv', np.float64, ('nv')) + + msftyyz_out = f_out.createVariable('msftyyz', np.float32, ('time', 'basin', 'z_l', 'yh'), fill_value=1.e20) + msftyyz_out.missing_value = 1.e20 + + msftyzsmpa_out = f_out.createVariable('msftyzsmpa', np.float32, ('time', 'basin', 'z_l', 'yh'), fill_value=1.e20) + msftyzsmpa_out.missing_value = 1.e20 + + msftyzmpa_out = f_out.createVariable('msftyzmpa', np.float32, ('time', 'basin', 'z_l', 'yh'), fill_value=1.e20) + msftyzmpa_out.missing_value = 1.e20 + + average_T1_out = f_out.createVariable('average_T1', np.float64, ('time')) + average_T2_out = f_out.createVariable('average_T2', np.float64, ('time')) + average_DT_out = f_out.createVariable('average_DT', np.float64, ('time')) + time_bnds_out = f_out.createVariable('time_bnds', np.float64, ('time', 'nv')) + + time_out.setncatts(tax.__dict__) + yh_out.setncatts(yh.__dict__) + z_l_out.setncatts(z_l.__dict__) + z_i_out.setncatts(z_i.__dict__) + nv_out.setncatts(nv.__dict__) + + for k in msftyyz.__dict__.keys(): + if k[0] != '_': msftyyz_out.setncattr(k,msftyyz.__dict__[k]) + + for k in msftyzsmpa.__dict__.keys(): + if k[0] != '_': msftyzsmpa_out.setncattr(k,msftyzsmpa.__dict__[k]) + + for k in msftyzmpa.__dict__.keys(): + if k[0] != '_': msftyzmpa_out.setncattr(k,msftyzmpa.__dict__[k]) + + average_T1_out.setncatts(average_T1.__dict__) + average_T2_out.setncatts(average_T2.__dict__) + average_DT_out.setncatts(average_DT.__dict__) + time_bnds_out.setncatts(time_bnds.__dict__) + + time_out[:] = np.array(tax[:]) + yh_out[:] = np.array(yh[:]) + z_l_out[:] = np.array(z_l[:]) + z_i_out[:] = np.array(z_i[:]) + nv_out[:] = np.array(nv[:]) + + msftyyz_out[:] = np.ma.array(msftyyz[:]) + msftyzsmpa_out[:] = np.ma.array(msftyzsmpa[:]) + msftyzmpa_out[:] = np.ma.array(msftyzmpa[:]) + + average_T1_out[:] = average_T1[:] + average_T2_out[:] = average_T2[:] + average_DT_out[:] = average_DT[:] + time_bnds_out[:] = time_bnds[:] + + region_out[:] = nc.stringtochar(region) + + f_out.close() + + exit(0) + +if __name__ == '__main__': + run() diff --git a/tools/analysis/m6toolbox.py b/tools/analysis/m6toolbox.py index ba44ae2252..7d925e6a85 100644 --- a/tools/analysis/m6toolbox.py +++ b/tools/analysis/m6toolbox.py @@ -169,6 +169,34 @@ def maskFromDepth(depth, zCellTop): wet[depth>-zCellTop] = 1 return wet +def MOCpsi(vh, vmsk=None): + """Sums 'vh' zonally and cumulatively in the vertical to yield an overturning stream function, psi(y,z).""" + shape = list(vh.shape); shape[-3] += 1 + psi = np.zeros(shape[:-1]) + if len(shape)==3: + for k in range(shape[-3]-1,0,-1): + if vmsk is None: psi[k-1,:] = psi[k,:] - vh[k-1].sum(axis=-1) + else: psi[k-1,:] = psi[k,:] - (vmsk*vh[k-1]).sum(axis=-1) + else: + for n in range(shape[0]): + for k in range(shape[-3]-1,0,-1): + if vmsk is None: psi[n,k-1,:] = psi[n,k,:] - vh[n,k-1].sum(axis=-1) + else: psi[n,k-1,:] = psi[n,k,:] - (vmsk*vh[n,k-1]).sum(axis=-1) + return psi + + +def moc_maskedarray(vh,mask=None): + if mask is not None: + _mask = np.ma.masked_where(np.not_equal(mask,1.),mask) + else: + _mask = 1. + _vh = np.ma.masked_where(np.ma.equal(vh,0.),vh) + _vh = _vh * _mask + _vh = _vh[:,::-1,:] # flip z-axis so running sum is from ocean floor to surface + _vh = np.ma.cumsum(np.ma.sum(_vh,axis=-1),axis=1) + _vh = _vh[:,::-1,:] # flip z-axis back to original order + return _vh + def nearestJI(x, y, xy0): """ Find (j,i) of cell with center nearest to (x0,y0). From 7173edbf2fc7f7629efa94a74ef1b6f237bdb997 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Wed, 20 Sep 2017 10:39:53 -0400 Subject: [PATCH 02/22] MOC reported on interfaces, not levels - z_i now used instead of z_l - MOC at the bottom-most interface is 0. by definition - topo masking appears to be correct --- tools/analysis/CMIP6_offline_diags.py | 18 +++++++++--------- tools/analysis/m6toolbox.py | 12 +++++++----- 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/tools/analysis/CMIP6_offline_diags.py b/tools/analysis/CMIP6_offline_diags.py index 5039c00cce..c22db96398 100755 --- a/tools/analysis/CMIP6_offline_diags.py +++ b/tools/analysis/CMIP6_offline_diags.py @@ -87,7 +87,7 @@ def main(args): #-- msftyyz varname = 'vmo' - msftyyz = np.ma.ones((len(tax),3,len(z_l),len(yh)))*0. + msftyyz = np.ma.ones((len(tax),3,len(z_i),len(yh)))*0. msftyyz[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) msftyyz[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) msftyyz[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) @@ -95,13 +95,13 @@ def main(args): msftyyz = np.ma.array(msftyyz,fill_value=1.e20) msftyyz.long_name = 'Ocean Y Overturning Mass Streamfunction' msftyyz.units = 'kg s-1' - msftyyz.cell_methods = 'z_l:sum yh:sum basin:mean time:mean' + msftyyz.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' msftyyz.time_avg_info = 'average_T1,average_T2,average_DT' msftyyz.standard_name = 'ocean_y_overturning_mass_streamfunction' #-- msftyzmpa varname = 'vhGM' - msftyzmpa = np.ma.ones((len(tax),3,len(z_l),len(yh)))*0. + msftyzmpa = np.ma.ones((len(tax),3,len(z_i),len(yh)))*0. msftyzmpa[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) msftyzmpa[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) msftyzmpa[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) @@ -109,14 +109,14 @@ def main(args): msftyzmpa = np.ma.array(msftyzmpa,fill_value=1.e20) msftyzmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized mesoscale advection' msftyzmpa.units = 'kg s-1' - msftyzmpa.cell_methods = 'z_l:sum yh:sum basin:mean time:mean' + msftyzmpa.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' msftyzmpa.time_avg_info = 'average_T1,average_T2,average_DT' msftyzmpa.standard_name = 'ocean_meridional_overturning_mass_streamfunction_due_to_parameterized_'+\ 'mesoscale_advection' #-- msftyzsmpa varname = 'vhml' - msftyzsmpa = np.ma.ones((len(tax),3,len(z_l),len(yh)))*0. + msftyzsmpa = np.ma.ones((len(tax),3,len(z_i),len(yh)))*0. msftyzsmpa[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) msftyzsmpa[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) msftyzsmpa[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) @@ -124,7 +124,7 @@ def main(args): msftyzsmpa = np.ma.array(msftyzsmpa,fill_value=1.e20) msftyzsmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized submesoscale advection' msftyzsmpa.units = 'kg s-1' - msftyzsmpa.cell_methods = 'z_l:sum yh:sum basin:mean time:mean' + msftyzsmpa.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' msftyzsmpa.time_avg_info = 'average_T1,average_T2,average_DT' msftyzsmpa.standard_name = 'ocean_meridional_overturning_mass_streamfunction_due_to_parameterized_'+\ 'submesoscale_advection' @@ -176,13 +176,13 @@ def main(args): z_i_out = f_out.createVariable('z_i', np.float64, ('z_i')) nv_out = f_out.createVariable('nv', np.float64, ('nv')) - msftyyz_out = f_out.createVariable('msftyyz', np.float32, ('time', 'basin', 'z_l', 'yh'), fill_value=1.e20) + msftyyz_out = f_out.createVariable('msftyyz', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=1.e20) msftyyz_out.missing_value = 1.e20 - msftyzsmpa_out = f_out.createVariable('msftyzsmpa', np.float32, ('time', 'basin', 'z_l', 'yh'), fill_value=1.e20) + msftyzsmpa_out = f_out.createVariable('msftyzsmpa', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=1.e20) msftyzsmpa_out.missing_value = 1.e20 - msftyzmpa_out = f_out.createVariable('msftyzmpa', np.float32, ('time', 'basin', 'z_l', 'yh'), fill_value=1.e20) + msftyzmpa_out = f_out.createVariable('msftyzmpa', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=1.e20) msftyzmpa_out.missing_value = 1.e20 average_T1_out = f_out.createVariable('average_T1', np.float64, ('time')) diff --git a/tools/analysis/m6toolbox.py b/tools/analysis/m6toolbox.py index 7d925e6a85..75e79870c4 100644 --- a/tools/analysis/m6toolbox.py +++ b/tools/analysis/m6toolbox.py @@ -190,11 +190,13 @@ def moc_maskedarray(vh,mask=None): _mask = np.ma.masked_where(np.not_equal(mask,1.),mask) else: _mask = 1. - _vh = np.ma.masked_where(np.ma.equal(vh,0.),vh) - _vh = _vh * _mask - _vh = _vh[:,::-1,:] # flip z-axis so running sum is from ocean floor to surface - _vh = np.ma.cumsum(np.ma.sum(_vh,axis=-1),axis=1) - _vh = _vh[:,::-1,:] # flip z-axis back to original order + _vh = vh * _mask + _vh_btm = np.ma.expand_dims(_vh[:,-1,:,:]*0.,axis=1) + _vh = np.ma.concatenate((_vh,_vh_btm),axis=1) + _vh = np.ma.sum(_vh,axis=-1) + _vh = _vh[:,::-1] # flip z-axis so running sum is from ocean floor to surface + _vh = np.ma.cumsum(_vh,axis=1) + _vh = _vh[:,::-1] # flip z-axis back to original order return _vh def nearestJI(x, y, xy0): From 52d52e3069ab52af780ceb99272c94a0d63ecd1a Mon Sep 17 00:00:00 2001 From: John Krasting Date: Tue, 26 Sep 2017 20:14:17 -0400 Subject: [PATCH 03/22] Consolidated refineDiag updates - Introduced "--refineDiagDir" command line option for saving output history file - Added refineDiag script for overturning diagnostics saved on rho vertical coordines (work-in-progress). - Added refineDiag script for heat transport diagnostics (work-in-progress) - Updated main MOM6_refineDiag.csh driver script --- tools/analysis/MOM6_refineDiag.csh | 4 + tools/analysis/m6toolbox.py | 2 +- tools/analysis/refineDiag_ocean_month.py | 218 ++++++++++++++++++ tools/analysis/refineDiag_ocean_month_rho2.py | 185 +++++++++++++++ ...e_diags.py => refineDiag_ocean_month_z.py} | 8 +- 5 files changed, 414 insertions(+), 3 deletions(-) create mode 100755 tools/analysis/refineDiag_ocean_month.py create mode 100755 tools/analysis/refineDiag_ocean_month_rho2.py rename tools/analysis/{CMIP6_offline_diags.py => refineDiag_ocean_month_z.py} (95%) diff --git a/tools/analysis/MOM6_refineDiag.csh b/tools/analysis/MOM6_refineDiag.csh index 25d2794493..ef3aa4fa79 100755 --- a/tools/analysis/MOM6_refineDiag.csh +++ b/tools/analysis/MOM6_refineDiag.csh @@ -100,6 +100,10 @@ if ( -e $yr1.ocean_static_no_mask_table.nc ) set ocean_static_file = $yr1.ocean_ $script_dir/EddyKineticEnergy.py -g $ocean_static_file -o $out_dir/refineDiag_ocean_annual/EddyKineticEnergy/EKE_mean_${yr1}.png -l ${yr1} $yr1.ocean_daily.nc $script_dir/calc_variance.py zos $yr1.ocean_daily.nc $refineDiagDir/$yr1.ocean_month_refined.nc +echo '==== Offline Diagnostics ====' +$script_dir/offline_overturning_z.py -b $ocean_static_file -r refineDiagDir $yr1.ocean_month_z.nc +$script_dir/offline_overturning_rho.py -b $ocean_static_file -r refineDiagDir $yr1.ocean_month_rho2.nc + echo " ---------- end yearly analysis ---------- " echo " -- end MOM6_refineDiag.csh -- " diff --git a/tools/analysis/m6toolbox.py b/tools/analysis/m6toolbox.py index 75e79870c4..eab1f7d5ea 100644 --- a/tools/analysis/m6toolbox.py +++ b/tools/analysis/m6toolbox.py @@ -193,7 +193,7 @@ def moc_maskedarray(vh,mask=None): _vh = vh * _mask _vh_btm = np.ma.expand_dims(_vh[:,-1,:,:]*0.,axis=1) _vh = np.ma.concatenate((_vh,_vh_btm),axis=1) - _vh = np.ma.sum(_vh,axis=-1) + _vh = np.ma.sum(_vh,axis=-1) * -1. _vh = _vh[:,::-1] # flip z-axis so running sum is from ocean floor to surface _vh = np.ma.cumsum(_vh,axis=1) _vh = _vh[:,::-1] # flip z-axis back to original order diff --git a/tools/analysis/refineDiag_ocean_month.py b/tools/analysis/refineDiag_ocean_month.py new file mode 100755 index 0000000000..238b52a01c --- /dev/null +++ b/tools/analysis/refineDiag_ocean_month.py @@ -0,0 +1,218 @@ +#!/usr/bin/env python + +import argparse +import m6toolbox +import netCDF4 as nc +import numpy as np +import os +import sys +import matplotlib.pyplot as plt + + + +##-- RefineDiag Script for CMIP6 +## +## Variables we intend to provide in z-coordinates: +## +## msftyyz -> vmo (ocean_z) * both 0.25 and 0.5 resolutions +## msftyzmpa -> vhGM (ocean_z) * applies only to 0.5 resolution +## msftyzsmpa -> vhml (ocean_z) * both 0.25 and 0.5 resolutions +## +## +## Variables we intend to provide in rho-coordinates: +## (potenital density referenced to 2000 m) +## +## msftyrho -> vmo +## msftyrhompa -> vhGM * applies only to 0.5 resolution +## +## +## 2-D variables we intent to provide: +## +## hfy -> T_ady_2d + ndiff_tracer_trans_y_2d_T * T_ady_2d needs to be converted to Watts (Cp = 3992.) +## ndiff_tracer_trans_y_2d_T already in Watts +## advective term in both 0.25 and 0.5 resolutions +## diffusive term only in 0.5 resolution +## +## hfx -> same recipie as above, expect for x-dimension +## hfbasin -> summed line of hfy in each basin +## +## +## Outstanding issues +## 1.) regirdding of vh, vhGM to rho-corrdinates +## 2.) vhGM and vhML units need to be in kg s-1 +## 2.) save out RHO_0 and Cp somewhere in netCDF files to key off of +## +## +## CMIP variables that will NOT be provided: +## +## hfbasinpmadv, hfbasinpsmadv, hfbasinpmdiff, hfbasinpadv +## (We advect the tracer with the residual mean velocity; individual terms cannot be diagnosed) +## +## htovgyre, htovovrt, sltovgyre, sltovovrt +## (Code for offline calculation not ready.) +## +##-- + +def run(): + parser = argparse.ArgumentParser(description='''CMIP6 RefineDiag Script for OM4''') + parser.add_argument('infile', type=str, help='''Input file''') + parser.add_argument('-b','--basinfile', type=str, default='', required=True, help='''File containing OM4 basin masks''') + parser.add_argument('-o','--outfile', type=str, default=None, help='''Output file name''') + parser.add_argument('-r','--refineDiagDir', type=str, default=None, help='''Path to refineDiagDir defined by FRE workflow)''') + args = parser.parse_args() + main(args) + +def main(args): + #-- Define Regions and their associated masks + # Note: The Atlantic should include other smaller bays/seas that are + # included in the definition used in meridional_overturning.py + + region = np.array(['atlantic_arctic_ocean','indian_pacific_ocean','global_ocean']) + + #-- Read basin masks + f_basin = nc.Dataset(args.basinfile) + basin_code = f_basin.variables['basin'][:] + + atlantic_arctic_mask = basin_code * 0. + atlantic_arctic_mask[(basin_code==2) | (basin_code==4) | (basin_code==6) | (basin_code==7) | (basin_code==8)] = 1. + + indo_pacific_mask = basin_code * 0. + indo_pacific_mask[(basin_code==3) | (basin_code==5)] = 1. + + #-- Read model data + f_in = nc.Dataset(args.infile) + + #-- Read in existing dimensions from history netcdf file + yq = f_in.variables['yq'] + xh = f_in.variables['xh'] + tax = f_in.variables['time'] + + #-- hfy + advective = f_in.variables['T_ady_2d'][:] + if 'ndiff_tracer_trans_y_2d_T' in f_in.variables.keys(): + diffusive = f_in.variables['ndiff_tracer_trans_y_2d_T'][:] + else: + print("Warning: diffusive term 'ndiff_tracer_trans_y_2d_T' not found. Check if this experiment is running with neutral diffusion.") + diffusive = advective * 0. + hfy = advective + diffusive + #hfy[hfy.mask] = 1.e20 + #hfy = np.ma.array(hfy,fill_value=1.e20) + hfy.long_name = 'Ocean Heat Y Transport' + hfy.units = 'W' + hfy.cell_methods = 'yq:point xh:mean time:mean' + hfy.time_avg_info = 'average_T1,average_T2,average_DT' + hfy.standard_name = 'ocean_heat_y_transport' + + #-- hfx + advective = f_in.variables['T_adx_2d'][:] + if 'ndiff_tracer_trans_x_2d_T' in f_in.variables.keys(): + diffusive = f_in.variables['ndiff_tracer_trans_x_2d_T'][:] + else: + print("Warning: diffusive term 'ndiff_tracer_trans_x_2d_T' not found. Check if this experiment is running with neutral diffusion.") + diffusive = advective * 0. + hfx = advective + diffusive + #hfx[hfx.mask] = 1.e20 + #hfx = np.ma.array(hfx,fill_value=1.e20) + hfx.long_name = 'Ocean Heat X Transport' + hfx.units = 'W' + hfx.cell_methods = 'yh:mean xq:point time:mean' + hfx.time_avg_info = 'average_T1,average_T2,average_DT' + hfx.standard_name = 'ocean_heat_x_transport' + + + #-- Read time bounds + nv = f_in.variables['nv'] + average_T1 = f_in.variables['average_T1'] + average_T2 = f_in.variables['average_T2'] + average_DT = f_in.variables['average_DT'] + time_bnds = f_in.variables['time_bnds'] + + #-- Generate output filename + if args.outfile is None: + if hasattr(f_in,'filename'): + args.outfile = f_in.filename + else: + args.outfile = os.path.basename(args.infile) + args.outfile = args.outfile.split('.') + args.outfile[-2] = args.outfile[-2]+'_refined' + args.outfile = '.'.join(args.outfile) + + if args.refineDiagDir is not None: + args.outfile = args.refineDiagDir + '/' + args.outfile + + #-- Write output file + try: + os.remove(args.outfile) + except: + pass + + if os.path.exists(args.outfile): + raise IOError('Output netCDF file already exists.') + exit(1) + + f_out = nc.Dataset(args.outfile, 'w', format='NETCDF3_CLASSIC') + f_out.setncatts(f_in.__dict__) + f_out.filename = args.outfile + + time_dim = f_out.createDimension('time', size=None) + basin_dim = f_out.createDimension('basin', size=3) + strlen_dim = f_out.createDimension('strlen', size=21) + yq_dim = f_out.createDimension('yq', size=len(yq[:])) + xh_dim = f_out.createDimension('xh', size=len(xh[:])) + nv_dim = f_out.createDimension('nv', size=len(nv[:])) + + time_out = f_out.createVariable('time', np.float64, ('time')) + #basin_out = f_out.createVariable('basin', np.int32, ('basin')) + yq_out = f_out.createVariable('yq', np.float64, ('yq')) + region_out = f_out.createVariable('region', 'c', ('basin', 'strlen')) + xh_out = f_out.createVariable('xh', np.float64, ('xh')) + nv_out = f_out.createVariable('nv', np.float64, ('nv')) + + hfy_out = f_out.createVariable('hfy', np.float32, ('time', 'yq', 'xh'), fill_value=1.e20) + hfy_out.missing_value = 1.e20 + + hfx_out = f_out.createVariable('hfx', np.float32, ('time', 'yq', 'xh'), fill_value=1.e20) + hfx_out.missing_value = 1.e20 + + average_T1_out = f_out.createVariable('average_T1', np.float64, ('time')) + average_T2_out = f_out.createVariable('average_T2', np.float64, ('time')) + average_DT_out = f_out.createVariable('average_DT', np.float64, ('time')) + time_bnds_out = f_out.createVariable('time_bnds', np.float64, ('time', 'nv')) + + time_out.setncatts(tax.__dict__) + yq_out.setncatts(yq.__dict__) + xh_out.setncatts(xh.__dict__) + nv_out.setncatts(nv.__dict__) + + for k in hfy.__dict__.keys(): + if k[0] != '_': hfy_out.setncattr(k,hfy.__dict__[k]) + + for k in hfx.__dict__.keys(): + if k[0] != '_': hfx_out.setncattr(k,hfx.__dict__[k]) + + average_T1_out.setncatts(average_T1.__dict__) + average_T2_out.setncatts(average_T2.__dict__) + average_DT_out.setncatts(average_DT.__dict__) + time_bnds_out.setncatts(time_bnds.__dict__) + + time_out[:] = np.array(tax[:]) + yq_out[:] = np.array(yq[:]) + xh_out[:] = np.array(xh[:]) + nv_out[:] = np.array(nv[:]) + + hfy_out[:] = np.ma.array(hfy[:]) + hfx_out[:] = np.ma.array(hfx[:]) + + average_T1_out[:] = average_T1[:] + average_T2_out[:] = average_T2[:] + average_DT_out[:] = average_DT[:] + time_bnds_out[:] = time_bnds[:] + + region_out[:] = nc.stringtochar(region) + + f_out.close() + + exit(0) + +if __name__ == '__main__': + run() diff --git a/tools/analysis/refineDiag_ocean_month_rho2.py b/tools/analysis/refineDiag_ocean_month_rho2.py new file mode 100755 index 0000000000..fa5e99c398 --- /dev/null +++ b/tools/analysis/refineDiag_ocean_month_rho2.py @@ -0,0 +1,185 @@ +#!/usr/bin/env python + +import argparse +import m6toolbox +import netCDF4 as nc +import numpy as np +import os +import sys +import matplotlib.pyplot as plt + +##-- RefineDiag Script for CMIP6 +## +## +## Variables we intend to provide in rho-coordinates: +## (potenital density referenced to 2000 m) +## +## msftyrho -> vmo +## msftyrhompa -> vhGM * applies only to 0.5 resolution +## +##-- + +def run(): + parser = argparse.ArgumentParser(description='''CMIP6 RefineDiag Script for OM4''') + parser.add_argument('infile', type=str, help='''Input file''') + parser.add_argument('-b','--basinfile', type=str, default='', required=True, help='''File containing OM4 basin masks''') + parser.add_argument('-o','--outfile', type=str, default=None, help='''Output file name''') + parser.add_argument('-r','--refineDiagDir', type=str, default=None, help='''Path to refineDiagDir defined by FRE workflow)''') + args = parser.parse_args() + main(args) + +def main(args): + #-- Define Regions and their associated masks + # Note: The Atlantic should include other smaller bays/seas that are + # included in the definition used in meridional_overturning.py + + region = np.array(['atlantic_arctic_ocean','indian_pacific_ocean','global_ocean']) + + #-- Read basin masks + f_basin = nc.Dataset(args.basinfile) + basin_code = f_basin.variables['basin'][:] + + atlantic_arctic_mask = basin_code * 0. + atlantic_arctic_mask[(basin_code==2) | (basin_code==4) | (basin_code==6) | (basin_code==7) | (basin_code==8)] = 1. + + indo_pacific_mask = basin_code * 0. + indo_pacific_mask[(basin_code==3) | (basin_code==5)] = 1. + + #-- Read model data + f_in = nc.Dataset(args.infile) + + #-- Read in existing dimensions from history netcdf file + yh = f_in.variables['yh'] + z_l = f_in.variables['z_l'] + z_i = f_in.variables['z_i'] + tax = f_in.variables['time'] + + #-- msftyrho + varname = 'vmo' + msftyrho = np.ma.ones((len(tax),3,len(z_i),len(yh)))*0. + msftyrho[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) + msftyrho[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) + msftyrho[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) + msftyrho[msftyrho.mask] = 1.e20 + msftyrho = np.ma.array(msftyrho,fill_value=1.e20) + msftyrho.long_name = 'Ocean Y Overturning Mass Streamfunction' + msftyrho.units = 'kg s-1' + msftyrho.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' + msftyrho.time_avg_info = 'average_T1,average_T2,average_DT' + msftyrho.standard_name = 'ocean_y_overturning_mass_streamfunction' + + #-- msftyrhompa + varname = 'vhGM' + msftyrhompa = np.ma.ones((len(tax),3,len(z_i),len(yh)))*0. + msftyrhompa[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) + msftyrhompa[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) + msftyrhompa[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) + msftyrhompa[msftyrhompa.mask] = 1.e20 + msftyrhompa = np.ma.array(msftyrhompa,fill_value=1.e20) + msftyrhompa.long_name = 'ocean Y overturning mass streamfunction due to parameterized mesoscale advection' + msftyrhompa.units = 'kg s-1' + msftyrhompa.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' + msftyrhompa.time_avg_info = 'average_T1,average_T2,average_DT' + msftyrhompa.standard_name = 'ocean_y_overturning_mass_streamfunction_due_to_parameterized_'+\ + 'mesoscale_advection' + + + #-- Read time bounds + nv = f_in.variables['nv'] + average_T1 = f_in.variables['average_T1'] + average_T2 = f_in.variables['average_T2'] + average_DT = f_in.variables['average_DT'] + time_bnds = f_in.variables['time_bnds'] + + #-- Generate output filename + if args.outfile is None: + if hasattr(f_in,'filename'): + args.outfile = f_in.filename + else: + args.outfile = os.path.basename(args.infile) + args.outfile = args.outfile.split('.') + args.outfile[-2] = args.outfile[-2]+'_refined' + args.outfile = '.'.join(args.outfile) + + if args.refineDiagDir is not None: + args.outfile = args.refineDiagDir + '/' + args.outfile + + #-- Write output file + try: + os.remove(args.outfile) + except: + pass + + if os.path.exists(args.outfile): + raise IOError('Output netCDF file already exists.') + exit(1) + + f_out = nc.Dataset(args.outfile, 'w', format='NETCDF3_CLASSIC') + f_out.setncatts(f_in.__dict__) + f_out.filename = args.outfile + + time_dim = f_out.createDimension('time', size=None) + basin_dim = f_out.createDimension('basin', size=3) + strlen_dim = f_out.createDimension('strlen', size=21) + yh_dim = f_out.createDimension('yh', size=len(yh[:])) + z_l_dim = f_out.createDimension('z_l', size=len(z_l[:])) + z_i_dim = f_out.createDimension('z_i', size=len(z_i[:])) + nv_dim = f_out.createDimension('nv', size=len(nv[:])) + + time_out = f_out.createVariable('time', np.float64, ('time')) + yh_out = f_out.createVariable('yh', np.float64, ('yh')) + region_out = f_out.createVariable('region', 'c', ('basin', 'strlen')) + z_l_out = f_out.createVariable('z_l', np.float64, ('z_l')) + z_i_out = f_out.createVariable('z_i', np.float64, ('z_i')) + nv_out = f_out.createVariable('nv', np.float64, ('nv')) + + msftyrho_out = f_out.createVariable('msftyrho', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=1.e20) + msftyrho_out.missing_value = 1.e20 + + msftyrhompa_out = f_out.createVariable('msftyrhompa', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=1.e20) + msftyrhompa_out.missing_value = 1.e20 + + average_T1_out = f_out.createVariable('average_T1', np.float64, ('time')) + average_T2_out = f_out.createVariable('average_T2', np.float64, ('time')) + average_DT_out = f_out.createVariable('average_DT', np.float64, ('time')) + time_bnds_out = f_out.createVariable('time_bnds', np.float64, ('time', 'nv')) + + time_out.setncatts(tax.__dict__) + yh_out.setncatts(yh.__dict__) + z_l_out.setncatts(z_l.__dict__) + z_i_out.setncatts(z_i.__dict__) + nv_out.setncatts(nv.__dict__) + + for k in msftyrho.__dict__.keys(): + if k[0] != '_': msftyrho_out.setncattr(k,msftyrho.__dict__[k]) + + for k in msftyrhompa.__dict__.keys(): + if k[0] != '_': msftyrhompa_out.setncattr(k,msftyrhompa.__dict__[k]) + + average_T1_out.setncatts(average_T1.__dict__) + average_T2_out.setncatts(average_T2.__dict__) + average_DT_out.setncatts(average_DT.__dict__) + time_bnds_out.setncatts(time_bnds.__dict__) + + time_out[:] = np.array(tax[:]) + yh_out[:] = np.array(yh[:]) + z_l_out[:] = np.array(z_l[:]) + z_i_out[:] = np.array(z_i[:]) + nv_out[:] = np.array(nv[:]) + + msftyrho_out[:] = np.ma.array(msftyrho[:]) + msftyrhompa_out[:] = np.ma.array(msftyrhompa[:]) + + average_T1_out[:] = average_T1[:] + average_T2_out[:] = average_T2[:] + average_DT_out[:] = average_DT[:] + time_bnds_out[:] = time_bnds[:] + + region_out[:] = nc.stringtochar(region) + + f_out.close() + + exit(0) + +if __name__ == '__main__': + run() diff --git a/tools/analysis/CMIP6_offline_diags.py b/tools/analysis/refineDiag_ocean_month_z.py similarity index 95% rename from tools/analysis/CMIP6_offline_diags.py rename to tools/analysis/refineDiag_ocean_month_z.py index c22db96398..54a674bb8d 100755 --- a/tools/analysis/CMIP6_offline_diags.py +++ b/tools/analysis/refineDiag_ocean_month_z.py @@ -55,7 +55,8 @@ def run(): parser = argparse.ArgumentParser(description='''CMIP6 RefineDiag Script for OM4''') parser.add_argument('infile', type=str, help='''Input file''') parser.add_argument('-b','--basinfile', type=str, default='', required=True, help='''File containing OM4 basin masks''') - parser.add_argument('-o','--outfile', type=str, default=None, help='''Output file name''') + parser.add_argument('-o','--outfile', type=str, default=None, help='''Output file name''') + parser.add_argument('-r','--refineDiagDir', type=str, default=None, help='''Path to refineDiagDir defined by FRE workflow)''') args = parser.parse_args() main(args) @@ -111,7 +112,7 @@ def main(args): msftyzmpa.units = 'kg s-1' msftyzmpa.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' msftyzmpa.time_avg_info = 'average_T1,average_T2,average_DT' - msftyzmpa.standard_name = 'ocean_meridional_overturning_mass_streamfunction_due_to_parameterized_'+\ + msftyzmpa.standard_name = 'ocean_y_overturning_mass_streamfunction_due_to_parameterized_'+\ 'mesoscale_advection' #-- msftyzsmpa @@ -146,6 +147,9 @@ def main(args): args.outfile[-2] = args.outfile[-2]+'_refined' args.outfile = '.'.join(args.outfile) + if args.refineDiagDir is not None: + args.outfile = args.refineDiagDir + '/' + args.outfile + #-- Write output file try: os.remove(args.outfile) From da8591d35bbf413cae5eee74935fd86ee565df35 Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Thu, 28 Sep 2017 18:27:11 -0400 Subject: [PATCH 04/22] Include wmo as part of refine diags A requested variable 'wmo' the upwards mass transport is a requested CMIP6 diagnostic which as of yet is not available from the prognostic model in z* space. Transport across a vertical interface in hybrid coordinates can be slightly different than what users may canonically think of as a upward (towards the surface) direction which is orthogonal to strictly meridional/zonal directions. To be consistent with this colloquial understanding of vertical direction, wmo is calculated from the convergence of umo and vmo which have been remapped to a z* coordinate. --- tools/analysis/refineDiag_ocean_month_z.py | 109 ++++++++++++--------- tools/analysis/verticalvelocity.py | 82 ++++++++++++++++ 2 files changed, 147 insertions(+), 44 deletions(-) create mode 100755 tools/analysis/verticalvelocity.py diff --git a/tools/analysis/refineDiag_ocean_month_z.py b/tools/analysis/refineDiag_ocean_month_z.py index 54a674bb8d..1b3e72f471 100755 --- a/tools/analysis/refineDiag_ocean_month_z.py +++ b/tools/analysis/refineDiag_ocean_month_z.py @@ -7,11 +7,12 @@ import os import sys import matplotlib.pyplot as plt +from verticalvelocity import calc_w_from_convergence ##-- RefineDiag Script for CMIP6 ## ## Variables we intend to provide in z-coordinates: -## +## ## msftyyz -> vmo (ocean_z) * both 0.25 and 0.5 resolutions ## msftyzmpa -> vhGM (ocean_z) * applies only to 0.5 resolution ## msftyzsmpa -> vhml (ocean_z) * both 0.25 and 0.5 resolutions @@ -20,7 +21,7 @@ ## Variables we intend to provide in rho-coordinates: ## (potenital density referenced to 2000 m) ## -## msftyrho -> vmo +## msftyrho -> vmo ## msftyrhompa -> vhGM * applies only to 0.5 resolution ## ## @@ -32,20 +33,20 @@ ## diffusive term only in 0.5 resolution ## ## hfx -> same recipie as above, expect for x-dimension -## hfbasin -> summed line of hfy in each basin -## +## hfbasin -> summed line of hfy in each basin +## ## ## Outstanding issues ## 1.) regirdding of vh, vhGM to rho-corrdinates ## 2.) vhGM and vhML units need to be in kg s-1 ## 2.) save out RHO_0 and Cp somewhere in netCDF files to key off of ## -## +## ## CMIP variables that will NOT be provided: ## ## hfbasinpmadv, hfbasinpsmadv, hfbasinpmdiff, hfbasinpadv ## (We advect the tracer with the residual mean velocity; individual terms cannot be diagnosed) -## +## ## htovgyre, htovovrt, sltovgyre, sltovovrt ## (Code for offline calculation not ready.) ## @@ -61,12 +62,13 @@ def run(): main(args) def main(args): + nc_misval = 1.e20 #-- Define Regions and their associated masks - # Note: The Atlantic should include other smaller bays/seas that are + # Note: The Atlantic should include other smaller bays/seas that are # included in the definition used in meridional_overturning.py - + region = np.array(['atlantic_arctic_ocean','indian_pacific_ocean','global_ocean']) - + #-- Read basin masks f_basin = nc.Dataset(args.basinfile) basin_code = f_basin.variables['basin'][:] @@ -76,24 +78,25 @@ def main(args): indo_pacific_mask = basin_code * 0. indo_pacific_mask[(basin_code==3) | (basin_code==5)] = 1. - + #-- Read model data f_in = nc.Dataset(args.infile) #-- Read in existing dimensions from history netcdf file + xh = f_in.variables['xh'] yh = f_in.variables['yh'] z_l = f_in.variables['z_l'] z_i = f_in.variables['z_i'] tax = f_in.variables['time'] - + #-- msftyyz varname = 'vmo' msftyyz = np.ma.ones((len(tax),3,len(z_i),len(yh)))*0. msftyyz[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) msftyyz[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) msftyyz[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) - msftyyz[msftyyz.mask] = 1.e20 - msftyyz = np.ma.array(msftyyz,fill_value=1.e20) + msftyyz[msftyyz.mask] = nc_misval + msftyyz = np.ma.array(msftyyz,fill_value=nc_misval) msftyyz.long_name = 'Ocean Y Overturning Mass Streamfunction' msftyyz.units = 'kg s-1' msftyyz.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' @@ -106,8 +109,8 @@ def main(args): msftyzmpa[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) msftyzmpa[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) msftyzmpa[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) - msftyzmpa[msftyzmpa.mask] = 1.e20 - msftyzmpa = np.ma.array(msftyzmpa,fill_value=1.e20) + msftyzmpa[msftyzmpa.mask] = nc_misval + msftyzmpa = np.ma.array(msftyzmpa,fill_value=nc_misval) msftyzmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized mesoscale advection' msftyzmpa.units = 'kg s-1' msftyzmpa.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' @@ -121,8 +124,8 @@ def main(args): msftyzsmpa[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) msftyzsmpa[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) msftyzsmpa[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) - msftyzsmpa[msftyzsmpa.mask] = 1.e20 - msftyzsmpa = np.ma.array(msftyzsmpa,fill_value=1.e20) + msftyzsmpa[msftyzsmpa.mask] = nc_misval + msftyzsmpa = np.ma.array(msftyzsmpa,fill_value=nc_misval) msftyzsmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized submesoscale advection' msftyzsmpa.units = 'kg s-1' msftyzsmpa.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' @@ -130,13 +133,24 @@ def main(args): msftyzsmpa.standard_name = 'ocean_meridional_overturning_mass_streamfunction_due_to_parameterized_'+\ 'submesoscale_advection' - #-- Read time bounds + #-- wmo + varname = 'wmo' + wmo = calc_w_from_convergence(infile, f_in.variables['umo'], f_in.variables['vmo']) + wmo[wmo.mask] = nc_misval + wmo = np.ma.array(wmo,fill_value=nc_misval) + wmo.long_name = 'Upward mass transport from resolved and parameterized advective transport' + wmo.units = 'kg s-1' + wmo.cell_methods = 'z_i:sum xh:sum yh:sum time:mean' + wmo.time_avg_info = 'average_T1,average_T2,average_DT' + wmo.standard_name = 'upward_ocean_mass_transport' + + #-- Read time bounds nv = f_in.variables['nv'] average_T1 = f_in.variables['average_T1'] average_T2 = f_in.variables['average_T2'] average_DT = f_in.variables['average_DT'] time_bnds = f_in.variables['time_bnds'] - + #-- Generate output filename if args.outfile is None: if hasattr(f_in,'filename'): @@ -146,7 +160,7 @@ def main(args): args.outfile = args.outfile.split('.') args.outfile[-2] = args.outfile[-2]+'_refined' args.outfile = '.'.join(args.outfile) - + if args.refineDiagDir is not None: args.outfile = args.refineDiagDir + '/' + args.outfile @@ -155,23 +169,24 @@ def main(args): os.remove(args.outfile) except: pass - + if os.path.exists(args.outfile): raise IOError('Output netCDF file already exists.') exit(1) - + f_out = nc.Dataset(args.outfile, 'w', format='NETCDF3_CLASSIC') f_out.setncatts(f_in.__dict__) f_out.filename = args.outfile - + time_dim = f_out.createDimension('time', size=None) basin_dim = f_out.createDimension('basin', size=3) strlen_dim = f_out.createDimension('strlen', size=21) + xh_dim = f_out.createDimension('xh', size=len(xh[:])) yh_dim = f_out.createDimension('yh', size=len(yh[:])) z_l_dim = f_out.createDimension('z_l', size=len(z_l[:])) z_i_dim = f_out.createDimension('z_i', size=len(z_i[:])) nv_dim = f_out.createDimension('nv', size=len(nv[:])) - + time_out = f_out.createVariable('time', np.float64, ('time')) #basin_out = f_out.createVariable('basin', np.int32, ('basin')) yh_out = f_out.createVariable('yh', np.float64, ('yh')) @@ -179,22 +194,26 @@ def main(args): z_l_out = f_out.createVariable('z_l', np.float64, ('z_l')) z_i_out = f_out.createVariable('z_i', np.float64, ('z_i')) nv_out = f_out.createVariable('nv', np.float64, ('nv')) - - msftyyz_out = f_out.createVariable('msftyyz', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=1.e20) - msftyyz_out.missing_value = 1.e20 - - msftyzsmpa_out = f_out.createVariable('msftyzsmpa', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=1.e20) - msftyzsmpa_out.missing_value = 1.e20 - - msftyzmpa_out = f_out.createVariable('msftyzmpa', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=1.e20) - msftyzmpa_out.missing_value = 1.e20 - + + msftyyz_out = f_out.createVariable('msftyyz', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=nc_misval) + msftyyz_out.missing_value = nc_misval + + msftyzsmpa_out = f_out.createVariable('msftyzsmpa', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=nc_misval) + msftyzsmpa_out.missing_value = nc_misval + + msftyzmpa_out = f_out.createVariable('msftyzmpa', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=nc_misval) + msftyzmpa_out.missing_value = nc_misval + + wmo_out = f_out.createVariable('wmo', np.float32, ('time', 'z_i', 'yh', 'xh'), fill_value=nc_misval) + wmo_out.missing_value = nc_misval + average_T1_out = f_out.createVariable('average_T1', np.float64, ('time')) average_T2_out = f_out.createVariable('average_T2', np.float64, ('time')) average_DT_out = f_out.createVariable('average_DT', np.float64, ('time')) time_bnds_out = f_out.createVariable('time_bnds', np.float64, ('time', 'nv')) - + time_out.setncatts(tax.__dict__) + xh_out.setncatts(xh.__dict__) yh_out.setncatts(yh.__dict__) z_l_out.setncatts(z_l.__dict__) z_i_out.setncatts(z_i.__dict__) @@ -202,37 +221,39 @@ def main(args): for k in msftyyz.__dict__.keys(): if k[0] != '_': msftyyz_out.setncattr(k,msftyyz.__dict__[k]) - + for k in msftyzsmpa.__dict__.keys(): if k[0] != '_': msftyzsmpa_out.setncattr(k,msftyzsmpa.__dict__[k]) - + for k in msftyzmpa.__dict__.keys(): if k[0] != '_': msftyzmpa_out.setncattr(k,msftyzmpa.__dict__[k]) - + average_T1_out.setncatts(average_T1.__dict__) average_T2_out.setncatts(average_T2.__dict__) average_DT_out.setncatts(average_DT.__dict__) time_bnds_out.setncatts(time_bnds.__dict__) - + time_out[:] = np.array(tax[:]) + xh_out[:] = np.array(xh[:]) yh_out[:] = np.array(yh[:]) z_l_out[:] = np.array(z_l[:]) z_i_out[:] = np.array(z_i[:]) nv_out[:] = np.array(nv[:]) - + msftyyz_out[:] = np.ma.array(msftyyz[:]) msftyzsmpa_out[:] = np.ma.array(msftyzsmpa[:]) msftyzmpa_out[:] = np.ma.array(msftyzmpa[:]) - + wmo_out[:] = np.ma.array(wmo[:]) + average_T1_out[:] = average_T1[:] average_T2_out[:] = average_T2[:] average_DT_out[:] = average_DT[:] time_bnds_out[:] = time_bnds[:] - + region_out[:] = nc.stringtochar(region) - + f_out.close() - + exit(0) if __name__ == '__main__': diff --git a/tools/analysis/verticalvelocity.py b/tools/analysis/verticalvelocity.py new file mode 100755 index 0000000000..b10f6e44ad --- /dev/null +++ b/tools/analysis/verticalvelocity.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python +# Estimate vertical mass transport (wmo) by the divergence of the horizontal mass transports. The vertical mass transport across an +# interface is the cumulative integral starting from the bottom of a water column. The sign convention is w>0 is an upward transport, +# (i.e., towards the surface of the ocean). By this convention, then div(u,v)<0 implies a negative (downward) transport and vice versa. +# Written by Andrew Shao (andrew.shao@noaa.gov) 29 September 2017 +import numpy as np + +def calc_w_from_convergence(infile, gridspec, u_var, v_var, wrapx = True, wrapy = False): + + tmax = u_var.shape[0] + + ntime, nk, nlat, nlon = u_var.shape + w = np.ma.zeros( (ntime, nk+1, nlat, nlon) ) + # Work timelevel by timelevel + for tidx in range(0,tmax): + # Get and process the u component + u_dat = u_var[tidx,:,:,:] + h_mask = np.logical_or(np.ma.getmask(u_dat), np.ma.getmask(np.roll(u_dat,1,axis=-1))) + u_dat = u_dat.filled(0.) + + # Get and process the v component + v_dat = v_var[tidx,:,:,:] + h_mask = np.logical_or(h_mask,np.ma.getmask(v_dat)) + h_mask = np.logical_or(h_mask,np.ma.getmask(np.roll(v_dat,1,axis=-2))) + v_dat = v_dat.filled(0.) + + # Order of subtraction based on upwind sign convention and desire for w>0 to correspond with upwards velocity + w[tidx,:-1,:,:] += u_dat-np.roll(u_dat,1,axis=-1) + if not wrapx: # If not wrapping, then convergence on westernmost side is simply so add back the rolled value + w[tidx,:-1,:,0] += u_dat[:-1,:,-1] + w[tidx,:-1,:,:] += v_dat-np.roll(v_dat,1,axis=-2) + if not wrapy: # If not wrapping, convergence on westernmost side is -v + w[tidx,:-1,0,:] += -v_dat[:,-1,:] + w[tidx,-1,:,:] = 0. + # Do a double-flip so that we integrate from the bottom + w[tidx,:-1,:,:] = w[tidx,-2::-1,:,:].cumsum(axis=0)[::-1,:,:] + # Mask if any of u[i-1], u[i], v[j-1], v[j] are not masked + w[tidx,:-1,:,:] = np.ma.masked_where(h_mask, w[tidx,:-1,:,:]) + # Bottom should always be zero, mask applied wherever the top interface is a valid value + w[tidx,-1,:,:] = np.ma.masked_where(h_mask[-2,:,:], w[tidx,-1,:,:]) + + return w + +if __name__ == "__main__": + try: import argparse + except: raise Exception('This version of python is not new enough. python 2.7 or newer is required.') + import netCDF4 + + parser = argparse.ArgumentParser(description= + '''Script for calculating vertical velocity from divergence of horizontal mass transports.''') + parser.add_argument('infile', type=str, help='''Daily file containing ssu,ssv.''') + parser.add_argument('--uname', type=str, default='umo', help='''Name of the u-component of mass transport''') + parser.add_argument('--vname', type=str, default='vmo', help='''Name of the v-component of mass transport''') + parser.add_argument('--wrapx', type=bool, default=True, help='''True if the x-component is reentrant''') + parser.add_argument('--wrapy', type=bool, default=False, help='''True if the y-component is reentrant''') + parser.add_argument('--gridspec', type=str, default ='', + help='''File containing variables wet,areacello. Usually the ocean_static.nc from diag_table.''') + parser.add_argument('--plot', type=bool, default=False, help='''Plot w at the 8th interface (kind of random)''') + cmdLineArgs = parser.parse_args() + + # Check for presence of variables + rootGroup = netCDF4.Dataset( cmdLineArgs.infile ) + if cmdLineArgs.uname not in rootGroup.variables: + raise Exception('Could not find "%s" in file "%s"' % (cmdLineArgs.uname, cmdLineArgs.infile)) + if cmdLineArgs.vname not in rootGroup.variables: + raise Exception('Could not find "%s" in file "%s"' % (cmdLineArgs.vname, cmdLineArgs.infile)) + + u_var = netCDF4.Dataset(cmdLineArgs.infile).variables[cmdLineArgs.uname] + v_var = netCDF4.Dataset(cmdLineArgs.infile).variables[cmdLineArgs.vname] + + w = calc_w_from_convergence( cmdLineArgs.infile, cmdLineArgs.gridspec, u_var, v_var, + cmdLineArgs.wrapx, cmdLineArgs.wrapy ) + if len(cmdLineArgs.gridspec) > 0: + area = netCDF4.Dataset(cmdLineArgs.gridspec).variables['areacello'][:,:] + w = w/(area*1035.0) + # Plot if requested + if cmdLineArgs.plot: + import matplotlib.pyplot as plt + plt.pcolormesh(w[0,7,:,:],vmin=-2e-5,vmax=2e-5,cmap=plt.cm.coolwarm) + plt.colorbar() + plt.show() + From a38e9834d035dc7e4cecf42f8e66d277ef8fda72 Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Thu, 28 Sep 2017 18:52:10 -0400 Subject: [PATCH 05/22] w had the wrong sign because divergence was calculated, not convergence --- tools/analysis/verticalvelocity.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tools/analysis/verticalvelocity.py b/tools/analysis/verticalvelocity.py index b10f6e44ad..c94cd0b873 100755 --- a/tools/analysis/verticalvelocity.py +++ b/tools/analysis/verticalvelocity.py @@ -25,11 +25,11 @@ def calc_w_from_convergence(infile, gridspec, u_var, v_var, wrapx = True, wrapy v_dat = v_dat.filled(0.) # Order of subtraction based on upwind sign convention and desire for w>0 to correspond with upwards velocity - w[tidx,:-1,:,:] += u_dat-np.roll(u_dat,1,axis=-1) - if not wrapx: # If not wrapping, then convergence on westernmost side is simply so add back the rolled value - w[tidx,:-1,:,0] += u_dat[:-1,:,-1] - w[tidx,:-1,:,:] += v_dat-np.roll(v_dat,1,axis=-2) - if not wrapy: # If not wrapping, convergence on westernmost side is -v + w[tidx,:-1,:,:] += np.roll(u_dat,1,axis=-1)-u_dat + if not wrapx: # If not wrapping, then convergence on westernmost side is simply so subtract back the rolled value + w[tidx,:-1,:,0] += -u_dat[:-1,:,-1] + w[tidx,:-1,:,:] += np.roll(v_dat,1,axis=-2)-v_dat + if not wrapy: # If not wrapping, convergence on westernmost side is v w[tidx,:-1,0,:] += -v_dat[:,-1,:] w[tidx,-1,:,:] = 0. # Do a double-flip so that we integrate from the bottom From 769eb6b208fd1ed3f10ffc2321b115b51304e031 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Sun, 1 Oct 2017 12:28:09 -0400 Subject: [PATCH 06/22] refineDiag script uses basin_codes.nc file - Also temporarily commented out daily variance of zos for testing purposes. --- tools/analysis/MOM6_refineDiag.csh | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/tools/analysis/MOM6_refineDiag.csh b/tools/analysis/MOM6_refineDiag.csh index ef3aa4fa79..b4b0ef1f38 100755 --- a/tools/analysis/MOM6_refineDiag.csh +++ b/tools/analysis/MOM6_refineDiag.csh @@ -91,18 +91,24 @@ ls -l set script_dir=${out_dir}/mom6/tools/analysis -echo '==Run some example annual scripts. These are not reviewed by scientists.' +set ocean_static_file = $yr1.ocean_static.nc +if ( -e $yr1.ocean_static_no_mask_table.nc ) set ocean_static_file = $yr1.ocean_static_no_mask_table.nc + +set basin_codes_file = $yr1.basin_codes.nc echo '====annual mean Eddy Kinetic Energy======' mkdir -p $out_dir/refineDiag_ocean_annual/EddyKineticEnergy -set ocean_static_file = $yr1.ocean_static.nc -if ( -e $yr1.ocean_static_no_mask_table.nc ) set ocean_static_file = $yr1.ocean_static_no_mask_table.nc $script_dir/EddyKineticEnergy.py -g $ocean_static_file -o $out_dir/refineDiag_ocean_annual/EddyKineticEnergy/EKE_mean_${yr1}.png -l ${yr1} $yr1.ocean_daily.nc -$script_dir/calc_variance.py zos $yr1.ocean_daily.nc $refineDiagDir/$yr1.ocean_month_refined.nc + +#-- ZOS temporarily commented out for now. This calculation needs to be merged into the +# refinedDiag_ocean_month.py script +#echo '====calculate zos as refineDiag====' +#$script_dir/calc_variance.py zos $yr1.ocean_daily.nc $refineDiagDir/$yr1.ocean_month_refined.nc echo '==== Offline Diagnostics ====' -$script_dir/offline_overturning_z.py -b $ocean_static_file -r refineDiagDir $yr1.ocean_month_z.nc -$script_dir/offline_overturning_rho.py -b $ocean_static_file -r refineDiagDir $yr1.ocean_month_rho2.nc +$script_dir/refineDiag_ocean_month.py -b $basin_codes_file -r refineDiagDir $yr1.ocean_month.nc +$script_dir/refineDiag_ocean_month_z.py -b $basin_codes_file -r refineDiagDir $yr1.ocean_month_z.nc +$script_dir/refineDiag_ocean_month_rho2.py -b $basin_codes_file -r refineDiagDir $yr1.ocean_month_rho2.nc echo " ---------- end yearly analysis ---------- " From 7f0ab79e346729780353548b6600b16f92c8a176 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Mon, 2 Oct 2017 13:32:02 -0400 Subject: [PATCH 07/22] Added hfbasin diagnostic to refineDiag script --- tools/analysis/refineDiag_ocean_month.py | 64 +++++++++++++----------- 1 file changed, 34 insertions(+), 30 deletions(-) diff --git a/tools/analysis/refineDiag_ocean_month.py b/tools/analysis/refineDiag_ocean_month.py index 238b52a01c..247001beb3 100755 --- a/tools/analysis/refineDiag_ocean_month.py +++ b/tools/analysis/refineDiag_ocean_month.py @@ -8,25 +8,9 @@ import sys import matplotlib.pyplot as plt - - ##-- RefineDiag Script for CMIP6 ## -## Variables we intend to provide in z-coordinates: -## -## msftyyz -> vmo (ocean_z) * both 0.25 and 0.5 resolutions -## msftyzmpa -> vhGM (ocean_z) * applies only to 0.5 resolution -## msftyzsmpa -> vhml (ocean_z) * both 0.25 and 0.5 resolutions -## -## -## Variables we intend to provide in rho-coordinates: -## (potenital density referenced to 2000 m) -## -## msftyrho -> vmo -## msftyrhompa -> vhGM * applies only to 0.5 resolution -## -## -## 2-D variables we intent to provide: +## 2-D variables we intend to provide: ## ## hfy -> T_ady_2d + ndiff_tracer_trans_y_2d_T * T_ady_2d needs to be converted to Watts (Cp = 3992.) ## ndiff_tracer_trans_y_2d_T already in Watts @@ -36,13 +20,6 @@ ## hfx -> same recipie as above, expect for x-dimension ## hfbasin -> summed line of hfy in each basin ## -## -## Outstanding issues -## 1.) regirdding of vh, vhGM to rho-corrdinates -## 2.) vhGM and vhML units need to be in kg s-1 -## 2.) save out RHO_0 and Cp somewhere in netCDF files to key off of -## -## ## CMIP variables that will NOT be provided: ## ## hfbasinpmadv, hfbasinpsmadv, hfbasinpmdiff, hfbasinpadv @@ -62,7 +39,20 @@ def run(): args = parser.parse_args() main(args) +def heat_trans_by_basin(x,mask=None): + if mask is not None: + varmask = np.sum(mask,axis=-1) + varmask = np.expand_dims(varmask,0) + varmask = np.where(np.equal(varmask,0.),True,False) + else: + mask = 1. + varmask = False + result = np.ma.sum((x*mask),axis=-1) + result.mask = varmask + return result + def main(args): + nc_misval = 1.e20 #-- Define Regions and their associated masks # Note: The Atlantic should include other smaller bays/seas that are # included in the definition used in meridional_overturning.py @@ -95,8 +85,6 @@ def main(args): print("Warning: diffusive term 'ndiff_tracer_trans_y_2d_T' not found. Check if this experiment is running with neutral diffusion.") diffusive = advective * 0. hfy = advective + diffusive - #hfy[hfy.mask] = 1.e20 - #hfy = np.ma.array(hfy,fill_value=1.e20) hfy.long_name = 'Ocean Heat Y Transport' hfy.units = 'W' hfy.cell_methods = 'yq:point xh:mean time:mean' @@ -111,14 +99,24 @@ def main(args): print("Warning: diffusive term 'ndiff_tracer_trans_x_2d_T' not found. Check if this experiment is running with neutral diffusion.") diffusive = advective * 0. hfx = advective + diffusive - #hfx[hfx.mask] = 1.e20 - #hfx = np.ma.array(hfx,fill_value=1.e20) hfx.long_name = 'Ocean Heat X Transport' hfx.units = 'W' hfx.cell_methods = 'yh:mean xq:point time:mean' hfx.time_avg_info = 'average_T1,average_T2,average_DT' hfx.standard_name = 'ocean_heat_x_transport' + #-- hfbasin + hfbasin = np.ma.ones((len(tax),3,len(yq)))*0. + hfbasin[:,0,:] = heat_trans_by_basin(hfy,mask=atlantic_arctic_mask) + hfbasin[:,1,:] = heat_trans_by_basin(hfy,mask=indo_pacific_mask) + hfbasin[:,2,:] = heat_trans_by_basin(hfy) + hfbasin[hfbasin.mask] = nc_misval + hfbasin = np.ma.array(hfbasin,fill_value=nc_misval) + hfbasin.long_name = 'Northward Ocean Heat Transport' + hfbasin.units = 'W' + hfbasin.cell_methods = 'xh:sum yq:sum basin:mean time:mean' + hfbasin.time_avg_info = 'average_T1,average_T2,average_DT' + hfbasin.standard_name = 'northward_ocean_heat_transport' #-- Read time bounds nv = f_in.variables['nv'] @@ -162,7 +160,6 @@ def main(args): nv_dim = f_out.createDimension('nv', size=len(nv[:])) time_out = f_out.createVariable('time', np.float64, ('time')) - #basin_out = f_out.createVariable('basin', np.int32, ('basin')) yq_out = f_out.createVariable('yq', np.float64, ('yq')) region_out = f_out.createVariable('region', 'c', ('basin', 'strlen')) xh_out = f_out.createVariable('xh', np.float64, ('xh')) @@ -174,6 +171,9 @@ def main(args): hfx_out = f_out.createVariable('hfx', np.float32, ('time', 'yq', 'xh'), fill_value=1.e20) hfx_out.missing_value = 1.e20 + hfbasin_out = f_out.createVariable('hfbasin', np.float32, ('time', 'basin', 'yq'), fill_value=1.e20) + hfbasin_out.missing_value = 1.e20 + average_T1_out = f_out.createVariable('average_T1', np.float64, ('time')) average_T2_out = f_out.createVariable('average_T2', np.float64, ('time')) average_DT_out = f_out.createVariable('average_DT', np.float64, ('time')) @@ -189,7 +189,10 @@ def main(args): for k in hfx.__dict__.keys(): if k[0] != '_': hfx_out.setncattr(k,hfx.__dict__[k]) - + + for k in hfbasin.__dict__.keys(): + if k[0] != '_': hfbasin_out.setncattr(k,hfbasin.__dict__[k]) + average_T1_out.setncatts(average_T1.__dict__) average_T2_out.setncatts(average_T2.__dict__) average_DT_out.setncatts(average_DT.__dict__) @@ -202,6 +205,7 @@ def main(args): hfy_out[:] = np.ma.array(hfy[:]) hfx_out[:] = np.ma.array(hfx[:]) + hfbasin_out[:] = np.ma.array(hfbasin[:]) average_T1_out[:] = average_T1[:] average_T2_out[:] = average_T2[:] From 5d039f628a9c78605a5855aefc84a81480d46aaf Mon Sep 17 00:00:00 2001 From: Andrew Shao Date: Tue, 3 Oct 2017 18:19:55 -0400 Subject: [PATCH 08/22] Adds a script to calculate transport in straits The diag table outputs diagnostic along some straits on a daily basis. The CMIP/OMIP requested 'mfo' diagnostic is a timeseries of monthly averaged, summed mass transport across the strait. The included script does this averaging and sums. The script can be extended if additional sections are added. --- tools/analysis/MOM6_refineDiag.csh | 35 ++--- tools/analysis/refineDiag_ocean_month_z.py | 24 +++- tools/analysis/strait_transport_timeseries.py | 124 ++++++++++++++++++ 3 files changed, 163 insertions(+), 20 deletions(-) create mode 100644 tools/analysis/strait_transport_timeseries.py diff --git a/tools/analysis/MOM6_refineDiag.csh b/tools/analysis/MOM6_refineDiag.csh index b4b0ef1f38..bde51f054a 100755 --- a/tools/analysis/MOM6_refineDiag.csh +++ b/tools/analysis/MOM6_refineDiag.csh @@ -2,23 +2,23 @@ #------------------------------------------------------------------------------ # MOM6_refineDiag.csh # -# DESCRIPTION: This is a script that is inteded to drive all the +# DESCRIPTION: This is a script that is inteded to drive all the # pre-postprocessing stages of data manipulations on analysis nodes. -# It is intended to be called by FRE at the "refineDiag" stage +# It is intended to be called by FRE at the "refineDiag" stage # which happens just before the components are post processed by frepp. -# To make this happen a path to the (would be) script should appear -# in the tag of the xmls, e.g., +# To make this happen a path to the (would be) script should appear +# in the tag of the xmls, e.g., # # Note that the above script should exist when frepp is called. -# This could be achieved by cloning the mom6 git repo in the section of the setup block -# in the corresponding the gfdl platfrom. E.g., +# This could be achieved by cloning the mom6 git repo in the section of the setup block +# in the corresponding the gfdl platfrom. E.g., # 0.: + z_i = v_vargroup['z_i'][:] + zidx = np.sum(z_i strait[sidx].zlim: + frac = min(1., (z_i[zidx] - strait[sidx].zlim)/(z_i[zidx]-z_i[zidx-1])) + vmo[:,-1,:] = vmo[:,-1,:]*frac + else: + vmo = v_vargroup['vmo'][:,:,:,:] + strait[sidx].transport += vmo.sum(axis=(1,2,3)) + Dataset(v_file[0]).close() + + if strait[sidx].is_meridional: + strait[sidx].time = u_vargroup['time'][:] + if strait[sidx].zlim > 0.: + z_i = u_vargroup['z_i'][:] + zidx = np.sum(z_i strait[sidx].zlim: + frac = min(1., (z_i[zidx] - strait[sidx].zlim)/(z_i[zidx]-z_i[zidx-1])) + umo[:,-1,:] = umo[:,-1,:]*frac + else: + umo = u_vargroup['umo'][:,:,:,:] + strait[sidx].transport += umo.sum(axis=(1,2,3)) + Dataset(u_file[0]).close() + if monthly_average: + strait[sidx].transport = make_monthly_averages(strait[sidx].transport) + strait[sidx].time = make_monthly_averages(strait[sidx].time) + ntime = strait[sidx].time.size + time = strait[sidx].time + + transport_array = np.zeros((ntime, nstraits)) + for sidx in range(0,nstraits): + transport_array[:,sidx] = strait[sidx].transport + return time, transport_array, strait + +def make_monthly_averages(data): + # This fundamentally assumes that the output is only one year long, any longer and we'd have to know the actual year to properly + # handle leap years + idx = daily_indices(data.size == 366) + idx_r = [range(0,idx[0])] ; + [idx_r.append(range(idx[t-1],idx[t])) for t in range(1,12)] + return np.array([np.mean(data[idx_r[t]]) for t in range(0,12)]) + +def daily_indices(leap = False): + days_in_months = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30 ,31, 30 ,31]) + if leap: # Change number of days in Feb if leapyear + days_in_months[1] = 29 + return days_in_months.cumsum() + + +def set_strait_info(): + + strait = [] + # The standard list of straits from CMIP/OMIP request + strait.append(strait_type("barents_opening", "Barents_opening", is_meridional = True)) + strait.append(strait_type("bering_strait", "Bering_Strait" , is_zonal = True)) + strait.append(strait_type("carribean_windward_passage", "Windward_Passage", is_zonal = True)) + strait.append(strait_type("davis_strait", "Davis_Strait", is_zonal = True)) + strait.append(strait_type("denmark_strait", "Denmark_Strait", is_zonal = True)) + strait.append(strait_type("drake_passage", "Drake_Passage", is_meridional = True)) + strait.append(strait_type("english_channel", "English_Channel", is_meridional = True)) + strait.append(strait_type("faroe_scotland_channel", "Faroe_Scotland", is_meridional = True)) + strait.append(strait_type("florida_bahamas_strait", "Florida_Bahamas", is_zonal = True)) + strait.append(strait_type("fram_strait", "Fram_Strait", is_zonal = True)) + strait.append(strait_type("gibraltar_strait", "Gibraltar_Strait", is_meridional = True)) + strait.append(strait_type("iceland_faroe_channel", "Iceland_Faroe", is_meridional = True, is_zonal = True)) + strait.append(strait_type("mozambique_channel", "Mozambique_Channel", is_zonal = True)) + strait.append(strait_type("pacific_equatorial_undercurrent", "Pacific_undercurrent", is_meridional = True, zlim = 350)) + strait.append(strait_type("taiwan_and_luzon_straits", "Taiwan_Luzon", is_meridional = True)) + strait.append(strait_type("agulhas_section", "Agulhas_section", is_meridional = True)) + strait.append(strait_type("indonesian_throughflow", "Indonesian_Throughflow", is_zonal = True)) + + return strait + +class strait_type: + def __init__(self, cmor_name, mom6_name, is_meridional = False, is_zonal = False, zlim = -1): + self.cmor_name = cmor_name + self.mom6_name = mom6_name + self.is_meridional = is_meridional + self.is_zonal = is_zonal + self.zlim = zlim + self.time = None + self.transport = None + + From 350462b8a6601315fe7166ba7005a3546e79a437 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Wed, 4 Oct 2017 08:53:21 -0400 Subject: [PATCH 09/22] Updated coordinates in rho2 refineDiag script --- tools/analysis/refineDiag_ocean_month_rho2.py | 42 +++++++++---------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/tools/analysis/refineDiag_ocean_month_rho2.py b/tools/analysis/refineDiag_ocean_month_rho2.py index fa5e99c398..f9a0122483 100755 --- a/tools/analysis/refineDiag_ocean_month_rho2.py +++ b/tools/analysis/refineDiag_ocean_month_rho2.py @@ -49,14 +49,14 @@ def main(args): f_in = nc.Dataset(args.infile) #-- Read in existing dimensions from history netcdf file - yh = f_in.variables['yh'] - z_l = f_in.variables['z_l'] - z_i = f_in.variables['z_i'] + yq = f_in.variables['yq'] + rho2_l = f_in.variables['rho2_l'] + rho2_i = f_in.variables['rho2_i'] tax = f_in.variables['time'] #-- msftyrho varname = 'vmo' - msftyrho = np.ma.ones((len(tax),3,len(z_i),len(yh)))*0. + msftyrho = np.ma.ones((len(tax),3,len(rho2_i),len(yq)))*0. msftyrho[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) msftyrho[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) msftyrho[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) @@ -64,13 +64,13 @@ def main(args): msftyrho = np.ma.array(msftyrho,fill_value=1.e20) msftyrho.long_name = 'Ocean Y Overturning Mass Streamfunction' msftyrho.units = 'kg s-1' - msftyrho.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' + msftyrho.cell_methods = 'rho2_i:sum yq:sum basin:mean time:mean' msftyrho.time_avg_info = 'average_T1,average_T2,average_DT' msftyrho.standard_name = 'ocean_y_overturning_mass_streamfunction' #-- msftyrhompa varname = 'vhGM' - msftyrhompa = np.ma.ones((len(tax),3,len(z_i),len(yh)))*0. + msftyrhompa = np.ma.ones((len(tax),3,len(rho2_i),len(yq)))*0. msftyrhompa[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) msftyrhompa[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) msftyrhompa[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) @@ -78,7 +78,7 @@ def main(args): msftyrhompa = np.ma.array(msftyrhompa,fill_value=1.e20) msftyrhompa.long_name = 'ocean Y overturning mass streamfunction due to parameterized mesoscale advection' msftyrhompa.units = 'kg s-1' - msftyrhompa.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' + msftyrhompa.cell_methods = 'rho2_i:sum yq:sum basin:mean time:mean' msftyrhompa.time_avg_info = 'average_T1,average_T2,average_DT' msftyrhompa.standard_name = 'ocean_y_overturning_mass_streamfunction_due_to_parameterized_'+\ 'mesoscale_advection' @@ -121,22 +121,22 @@ def main(args): time_dim = f_out.createDimension('time', size=None) basin_dim = f_out.createDimension('basin', size=3) strlen_dim = f_out.createDimension('strlen', size=21) - yh_dim = f_out.createDimension('yh', size=len(yh[:])) - z_l_dim = f_out.createDimension('z_l', size=len(z_l[:])) - z_i_dim = f_out.createDimension('z_i', size=len(z_i[:])) + yq_dim = f_out.createDimension('yq', size=len(yq[:])) + rho2_l_dim = f_out.createDimension('rho2_l', size=len(rho2_l[:])) + rho2_i_dim = f_out.createDimension('rho2_i', size=len(rho2_i[:])) nv_dim = f_out.createDimension('nv', size=len(nv[:])) time_out = f_out.createVariable('time', np.float64, ('time')) - yh_out = f_out.createVariable('yh', np.float64, ('yh')) + yq_out = f_out.createVariable('yq', np.float64, ('yq')) region_out = f_out.createVariable('region', 'c', ('basin', 'strlen')) - z_l_out = f_out.createVariable('z_l', np.float64, ('z_l')) - z_i_out = f_out.createVariable('z_i', np.float64, ('z_i')) + rho2_l_out = f_out.createVariable('rho2_l', np.float64, ('rho2_l')) + rho2_i_out = f_out.createVariable('rho2_i', np.float64, ('rho2_i')) nv_out = f_out.createVariable('nv', np.float64, ('nv')) - msftyrho_out = f_out.createVariable('msftyrho', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=1.e20) + msftyrho_out = f_out.createVariable('msftyrho', np.float32, ('time', 'basin', 'rho2_i', 'yq'), fill_value=1.e20) msftyrho_out.missing_value = 1.e20 - msftyrhompa_out = f_out.createVariable('msftyrhompa', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=1.e20) + msftyrhompa_out = f_out.createVariable('msftyrhompa', np.float32, ('time', 'basin', 'rho2_i', 'yq'), fill_value=1.e20) msftyrhompa_out.missing_value = 1.e20 average_T1_out = f_out.createVariable('average_T1', np.float64, ('time')) @@ -145,9 +145,9 @@ def main(args): time_bnds_out = f_out.createVariable('time_bnds', np.float64, ('time', 'nv')) time_out.setncatts(tax.__dict__) - yh_out.setncatts(yh.__dict__) - z_l_out.setncatts(z_l.__dict__) - z_i_out.setncatts(z_i.__dict__) + yq_out.setncatts(yq.__dict__) + rho2_l_out.setncatts(rho2_l.__dict__) + rho2_i_out.setncatts(rho2_i.__dict__) nv_out.setncatts(nv.__dict__) for k in msftyrho.__dict__.keys(): @@ -162,9 +162,9 @@ def main(args): time_bnds_out.setncatts(time_bnds.__dict__) time_out[:] = np.array(tax[:]) - yh_out[:] = np.array(yh[:]) - z_l_out[:] = np.array(z_l[:]) - z_i_out[:] = np.array(z_i[:]) + yq_out[:] = np.array(yq[:]) + rho2_l_out[:] = np.array(rho2_l[:]) + rho2_i_out[:] = np.array(rho2_i[:]) nv_out[:] = np.array(nv[:]) msftyrho_out[:] = np.ma.array(msftyrho[:]) From b235c772e955c0d156210e2b30ba37b2618b1eb0 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Wed, 4 Oct 2017 10:51:06 -0400 Subject: [PATCH 10/22] Masked Indo-Pacific heat transport below 34S - Southern latititude boudaries different between the Indian and Pacific basins. This shows up as a discontinuity in the Y-ward ocean heat transport when summing the two basins. - Data masked below 34 S and comment added to the 'hfbasin' diagnostic --- tools/analysis/refineDiag_ocean_month.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/tools/analysis/refineDiag_ocean_month.py b/tools/analysis/refineDiag_ocean_month.py index 247001beb3..4c5db85bb4 100755 --- a/tools/analysis/refineDiag_ocean_month.py +++ b/tools/analysis/refineDiag_ocean_month.py @@ -39,7 +39,7 @@ def run(): args = parser.parse_args() main(args) -def heat_trans_by_basin(x,mask=None): +def heat_trans_by_basin(x,mask=None,lat=None,minlat=None): if mask is not None: varmask = np.sum(mask,axis=-1) varmask = np.expand_dims(varmask,0) @@ -49,6 +49,10 @@ def heat_trans_by_basin(x,mask=None): varmask = False result = np.ma.sum((x*mask),axis=-1) result.mask = varmask + if (minlat is not None) and (lat is not None): + if len(lat.shape) == 1: + lat = np.expand_dims(lat,0) + result = np.ma.masked_where(np.less(lat,minlat),result) return result def main(args): @@ -68,7 +72,7 @@ def main(args): indo_pacific_mask = basin_code * 0. indo_pacific_mask[(basin_code==3) | (basin_code==5)] = 1. - + #-- Read model data f_in = nc.Dataset(args.infile) @@ -108,13 +112,14 @@ def main(args): #-- hfbasin hfbasin = np.ma.ones((len(tax),3,len(yq)))*0. hfbasin[:,0,:] = heat_trans_by_basin(hfy,mask=atlantic_arctic_mask) - hfbasin[:,1,:] = heat_trans_by_basin(hfy,mask=indo_pacific_mask) + hfbasin[:,1,:] = heat_trans_by_basin(hfy,mask=indo_pacific_mask,minlat=-34,lat=yq) hfbasin[:,2,:] = heat_trans_by_basin(hfy) hfbasin[hfbasin.mask] = nc_misval hfbasin = np.ma.array(hfbasin,fill_value=nc_misval) hfbasin.long_name = 'Northward Ocean Heat Transport' hfbasin.units = 'W' hfbasin.cell_methods = 'xh:sum yq:sum basin:mean time:mean' + hfbasin.comment = 'Indo-Pacific heat transport begins at 34 S' hfbasin.time_avg_info = 'average_T1,average_T2,average_DT' hfbasin.standard_name = 'northward_ocean_heat_transport' From 6ff7aac8790304ad33a6de4025699c8105aec59f Mon Sep 17 00:00:00 2001 From: John Krasting Date: Wed, 4 Oct 2017 10:56:27 -0400 Subject: [PATCH 11/22] Modified interface to calc_w_from_convergence() - infile and gridspec arguments appear to be unused. - temporarily commented out mfo diagnostic ; failing after merge --- tools/analysis/refineDiag_ocean_month_z.py | 70 +++++++++++++--------- tools/analysis/verticalvelocity.py | 4 +- 2 files changed, 44 insertions(+), 30 deletions(-) diff --git a/tools/analysis/refineDiag_ocean_month_z.py b/tools/analysis/refineDiag_ocean_month_z.py index cf54c3fb06..4d2c0b279a 100755 --- a/tools/analysis/refineDiag_ocean_month_z.py +++ b/tools/analysis/refineDiag_ocean_month_z.py @@ -86,13 +86,23 @@ def main(args): #-- Read in existing dimensions from history netcdf file xh = f_in.variables['xh'] yh = f_in.variables['yh'] + yq = f_in.variables['yq'] z_l = f_in.variables['z_l'] z_i = f_in.variables['z_i'] tax = f_in.variables['time'] + #-- Note: based on conversations with @adcroft, the overturning should be reported on the interfaces, z_i. + # Also, the nominal latitude is insufficient for the basin-average fields. Based on the methods in + # meridional_overturning.py, the latitude dimension should be: + # + # y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][::2,::2] + # yy = y[1:,:].max(axis=-1)+0*z + # + # The quanity 'yy' above is numerically-equivalent to 'yq' + #-- msftyyz varname = 'vmo' - msftyyz = np.ma.ones((len(tax),3,len(z_i),len(yh)))*0. + msftyyz = np.ma.ones((len(tax),3,len(z_i),len(yq)))*0. msftyyz[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) msftyyz[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) msftyyz[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) @@ -100,13 +110,13 @@ def main(args): msftyyz = np.ma.array(msftyyz,fill_value=nc_misval) msftyyz.long_name = 'Ocean Y Overturning Mass Streamfunction' msftyyz.units = 'kg s-1' - msftyyz.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' + msftyyz.cell_methods = 'z_i:sum yq:sum basin:mean time:mean' msftyyz.time_avg_info = 'average_T1,average_T2,average_DT' msftyyz.standard_name = 'ocean_y_overturning_mass_streamfunction' #-- msftyzmpa varname = 'vhGM' - msftyzmpa = np.ma.ones((len(tax),3,len(z_i),len(yh)))*0. + msftyzmpa = np.ma.ones((len(tax),3,len(z_i),len(yq)))*0. msftyzmpa[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) msftyzmpa[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) msftyzmpa[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) @@ -114,14 +124,14 @@ def main(args): msftyzmpa = np.ma.array(msftyzmpa,fill_value=nc_misval) msftyzmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized mesoscale advection' msftyzmpa.units = 'kg s-1' - msftyzmpa.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' + msftyzmpa.cell_methods = 'z_i:sum yq:sum basin:mean time:mean' msftyzmpa.time_avg_info = 'average_T1,average_T2,average_DT' msftyzmpa.standard_name = 'ocean_y_overturning_mass_streamfunction_due_to_parameterized_'+\ 'mesoscale_advection' #-- msftyzsmpa varname = 'vhml' - msftyzsmpa = np.ma.ones((len(tax),3,len(z_i),len(yh)))*0. + msftyzsmpa = np.ma.ones((len(tax),3,len(z_i),len(yq)))*0. msftyzsmpa[:,0,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=atlantic_arctic_mask) msftyzsmpa[:,1,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:],mask=indo_pacific_mask) msftyzsmpa[:,2,:,:] = m6toolbox.moc_maskedarray(f_in.variables[varname][:]) @@ -129,14 +139,14 @@ def main(args): msftyzsmpa = np.ma.array(msftyzsmpa,fill_value=nc_misval) msftyzsmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized submesoscale advection' msftyzsmpa.units = 'kg s-1' - msftyzsmpa.cell_methods = 'z_i:sum yh:sum basin:mean time:mean' + msftyzsmpa.cell_methods = 'z_i:sum yq:sum basin:mean time:mean' msftyzsmpa.time_avg_info = 'average_T1,average_T2,average_DT' msftyzsmpa.standard_name = 'ocean_meridional_overturning_mass_streamfunction_due_to_parameterized_'+\ 'submesoscale_advection' #-- wmo varname = 'wmo' - wmo = calc_w_from_convergence(infile, f_in.variables['umo'], f_in.variables['vmo']) + wmo = calc_w_from_convergence(f_in.variables['umo'], f_in.variables['vmo']) wmo[wmo.mask] = nc_misval wmo = np.ma.array(wmo,fill_value=nc_misval) wmo.long_name = 'Upward mass transport from resolved and parameterized advective transport' @@ -146,15 +156,15 @@ def main(args): wmo.standard_name = 'upward_ocean_mass_transport' #-- mfo - _, mfo, straits = sum_transport_in_straits(args.straitdir, monthly_average = True) - mfo[mfo].mask] = nc_misval - mfo = np.ma.array(mfo, fill_value=nc_misval) - strait_names = np.array( [strait.cmor_name for strait in straits] ) - mfo.long_name = 'Sea Water Transport' - mfo.units = 'kg s-1' - mfo.cell_methods = 'time:mean' - mfo.time_avg_info = 'average_T1,average_T2,average_DT' - mfo.standard_name = 'sea_water_transport_across_line' + #_, mfo, straits = sum_transport_in_straits(args.straitdir, monthly_average = True) + #mfo[mfo].mask = nc_misval + #mfo = np.ma.array(mfo, fill_value=nc_misval) + #strait_names = np.array( [strait.cmor_name for strait in straits] ) + #mfo.long_name = 'Sea Water Transport' + #mfo.units = 'kg s-1' + #mfo.cell_methods = 'time:mean' + #mfo.time_avg_info = 'average_T1,average_T2,average_DT' + #mfo.standard_name = 'sea_water_transport_across_line' #-- Read time bounds nv = f_in.variables['nv'] @@ -192,10 +202,11 @@ def main(args): time_dim = f_out.createDimension('time', size=None) basin_dim = f_out.createDimension('basin', size=3) - strait_dim = f_out.createDimension('strait', size=len(straits)) - strlen_dim = f_out.createDimension('strlen', size=31) + #strait_dim = f_out.createDimension('strait', size=len(straits)) + strlen_dim = f_out.createDimension('strlen', size=21) xh_dim = f_out.createDimension('xh', size=len(xh[:])) yh_dim = f_out.createDimension('yh', size=len(yh[:])) + yq_dim = f_out.createDimension('yq', size=len(yq[:])) z_l_dim = f_out.createDimension('z_l', size=len(z_l[:])) z_i_dim = f_out.createDimension('z_i', size=len(z_i[:])) nv_dim = f_out.createDimension('nv', size=len(nv[:])) @@ -203,26 +214,27 @@ def main(args): time_out = f_out.createVariable('time', np.float64, ('time')) #basin_out = f_out.createVariable('basin', np.int32, ('basin')) yh_out = f_out.createVariable('yh', np.float64, ('yh')) + yq_out = f_out.createVariable('yq', np.float64, ('yq')) region_out = f_out.createVariable('region', 'c', ('basin', 'strlen')) - strait_out = f_out.createVariable('strait', 'c', ('strait', 'strlen')) + #strait_out = f_out.createVariable('strait', 'c', ('strait', 'strlen')) z_l_out = f_out.createVariable('z_l', np.float64, ('z_l')) z_i_out = f_out.createVariable('z_i', np.float64, ('z_i')) nv_out = f_out.createVariable('nv', np.float64, ('nv')) - msftyyz_out = f_out.createVariable('msftyyz', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=nc_misval) + msftyyz_out = f_out.createVariable('msftyyz', np.float32, ('time', 'basin', 'z_i', 'yq'), fill_value=nc_misval) msftyyz_out.missing_value = nc_misval - msftyzsmpa_out = f_out.createVariable('msftyzsmpa', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=nc_misval) + msftyzsmpa_out = f_out.createVariable('msftyzsmpa', np.float32, ('time', 'basin', 'z_i', 'yq'), fill_value=nc_misval) msftyzsmpa_out.missing_value = nc_misval - msftyzmpa_out = f_out.createVariable('msftyzmpa', np.float32, ('time', 'basin', 'z_i', 'yh'), fill_value=nc_misval) + msftyzmpa_out = f_out.createVariable('msftyzmpa', np.float32, ('time', 'basin', 'z_i', 'yq'), fill_value=nc_misval) msftyzmpa_out.missing_value = nc_misval wmo_out = f_out.createVariable('wmo', np.float32, ('time', 'z_i', 'yh', 'xh'), fill_value=nc_misval) wmo_out.missing_value = nc_misval - mfo_out = f_out.createVariable('wmo', np.float32, ('time', 'strait'), fill_value=nc_misval) - mfo_out.missing_value = nc_misval + #mfo_out = f_out.createVariable('wmo', np.float32, ('time', 'strait'), fill_value=nc_misval) + #mfo_out.missing_value = nc_misval average_T1_out = f_out.createVariable('average_T1', np.float64, ('time')) average_T2_out = f_out.createVariable('average_T2', np.float64, ('time')) @@ -230,8 +242,9 @@ def main(args): time_bnds_out = f_out.createVariable('time_bnds', np.float64, ('time', 'nv')) time_out.setncatts(tax.__dict__) - xh_out.setncatts(xh.__dict__) + #xh_out.setncatts(xh.__dict__) yh_out.setncatts(yh.__dict__) + yq_out.setncatts(yq.__dict__) z_l_out.setncatts(z_l.__dict__) z_i_out.setncatts(z_i.__dict__) nv_out.setncatts(nv.__dict__) @@ -251,8 +264,9 @@ def main(args): time_bnds_out.setncatts(time_bnds.__dict__) time_out[:] = np.array(tax[:]) - xh_out[:] = np.array(xh[:]) + #xh_out[:] = np.array(xh[:]) yh_out[:] = np.array(yh[:]) + yq_out[:] = np.array(yq[:]) z_l_out[:] = np.array(z_l[:]) z_i_out[:] = np.array(z_i[:]) nv_out[:] = np.array(nv[:]) @@ -261,7 +275,7 @@ def main(args): msftyzsmpa_out[:] = np.ma.array(msftyzsmpa[:]) msftyzmpa_out[:] = np.ma.array(msftyzmpa[:]) wmo_out[:] = np.ma.array(wmo[:]) - mfo_out[:] = np.ma.array(mfo[:]) + #mfo_out[:] = np.ma.array(mfo[:]) average_T1_out[:] = average_T1[:] average_T2_out[:] = average_T2[:] @@ -269,7 +283,7 @@ def main(args): time_bnds_out[:] = time_bnds[:] region_out[:] = nc.stringtochar(region) - strait_out[:] = nc.stringtochar(strait_names) + #strait_out[:] = nc.stringtochar(strait_names) f_out.close() exit(0) diff --git a/tools/analysis/verticalvelocity.py b/tools/analysis/verticalvelocity.py index c94cd0b873..451c1c4c07 100755 --- a/tools/analysis/verticalvelocity.py +++ b/tools/analysis/verticalvelocity.py @@ -5,7 +5,7 @@ # Written by Andrew Shao (andrew.shao@noaa.gov) 29 September 2017 import numpy as np -def calc_w_from_convergence(infile, gridspec, u_var, v_var, wrapx = True, wrapy = False): +def calc_w_from_convergence(u_var, v_var, wrapx = True, wrapy = False): tmax = u_var.shape[0] @@ -68,7 +68,7 @@ def calc_w_from_convergence(infile, gridspec, u_var, v_var, wrapx = True, wrapy u_var = netCDF4.Dataset(cmdLineArgs.infile).variables[cmdLineArgs.uname] v_var = netCDF4.Dataset(cmdLineArgs.infile).variables[cmdLineArgs.vname] - w = calc_w_from_convergence( cmdLineArgs.infile, cmdLineArgs.gridspec, u_var, v_var, + w = calc_w_from_convergence( u_var, v_var, cmdLineArgs.wrapx, cmdLineArgs.wrapy ) if len(cmdLineArgs.gridspec) > 0: area = netCDF4.Dataset(cmdLineArgs.gridspec).variables['areacello'][:,:] From 69f35ec016de3554c678039362cc65d2307b912b Mon Sep 17 00:00:00 2001 From: John Krasting Date: Wed, 4 Oct 2017 11:40:51 -0400 Subject: [PATCH 12/22] Updated calc_variance.py - Script now either created the ocean_month_refined.nc file if it does not exist, or it appends variables to the file if does exist. - For CMIP6 runs, recommended that this script be called last in the refineDiag call order for best compatibility --- tools/analysis/MOM6_refineDiag.csh | 8 +-- tools/analysis/calc_variance.py | 102 ++++++++++++++++++++--------- 2 files changed, 76 insertions(+), 34 deletions(-) diff --git a/tools/analysis/MOM6_refineDiag.csh b/tools/analysis/MOM6_refineDiag.csh index bde51f054a..d1fe08d6c4 100755 --- a/tools/analysis/MOM6_refineDiag.csh +++ b/tools/analysis/MOM6_refineDiag.csh @@ -101,15 +101,15 @@ echo '====annual mean Eddy Kinetic Energy======' mkdir -p $out_dir/refineDiag_ocean_annual/EddyKineticEnergy $script_dir/EddyKineticEnergy.py -g $ocean_static_file -o $out_dir/refineDiag_ocean_annual/EddyKineticEnergy/EKE_mean_${yr1}.png -l ${yr1} $yr1.ocean_daily.nc -#-- ZOS temporarily commented out for now. This calculation needs to be merged into the -# refinedDiag_ocean_month.py script -#echo '====calculate zos as refineDiag====' -#$script_dir/calc_variance.py zos $yr1.ocean_daily.nc $refineDiagDir/$yr1.ocean_month_refined.nc echo '==== Offline Diagnostics ====' $script_dir/refineDiag_ocean_month.py -b $basin_codes_file -r refineDiagDir $yr1.ocean_month.nc $script_dir/refineDiag_ocean_month_z.py -b $basin_codes_file -r refineDiagDir -s $straitdir $yr1.ocean_month_z.nc $script_dir/refineDiag_ocean_month_rho2.py -b $basin_codes_file -r refineDiagDir $yr1.ocean_month_rho2.nc +$script_dir/calc_variance.py zos $yr1.ocean_daily.nc $refineDiagDir/$yr1.ocean_month_refined.nc + +#-- Note: the calc_variance script should be called LAST, since it appends to the ocean_month_refined.nc file +# that is created first by the other scripts. echo " ---------- end yearly analysis ---------- " diff --git a/tools/analysis/calc_variance.py b/tools/analysis/calc_variance.py index 1c4d84f46c..378fb6f890 100755 --- a/tools/analysis/calc_variance.py +++ b/tools/analysis/calc_variance.py @@ -2,6 +2,7 @@ import netCDF4 import numpy +import os try: import argparse except: raise Exception('This version of python is not new enough. python 2.7 or newer is required.') @@ -23,7 +24,25 @@ nxy = shape[1:] if args.verbose: print 'Creating',args.annual_file -nc_out = netCDF4.Dataset( args.annual_file, 'w', format='NETCDF3_64BIT' ) +if os.path.exists(args.annual_file): + mode = 'r+' + append = True +else: + mode = 'w' + append = False +nc_out = netCDF4.Dataset( args.annual_file, mode, format='NETCDF3_CLASSIC' ) + +if append is True: + if 'time' in nc_out.variables.keys(): + time_var = 'time' + elif 'Time' in nc_out.variables.keys(): + time_var = 'Time' + else: + nc_out.close() + raise ValueError('Time dimension could not be found in existing file.') + if len(nc_out.variables[time_var][:]) <= 1: + nc_out.close() + raise ValueError('Existing file has only one time value. Assuming this is annual output. Aborting.') if nt==365: days_in_month = [31,28,31,30,31,30,31,31,30,31,30,31] elif nt==366: days_in_month = [31,29,31,30,31,30,31,31,30,31,30,31] @@ -32,12 +51,15 @@ # Create dimensions for d in variable.dimensions: - if nc_in.dimensions[d].isunlimited(): nc_out.createDimension(d, None) - else: nc_out.createDimension(d, len(nc_in.dimensions[d])) + if d not in nc_out.dimensions: + if nc_in.dimensions[d].isunlimited(): nc_out.createDimension(d, None) + else: nc_out.createDimension(d, len(nc_in.dimensions[d])) # Copy global attributes -for a in nc_in.ncattrs(): - nc_out.__setattr__(a ,nc_in.__getattr__(a)) +if append is False: + for a in nc_in.ncattrs(): + if a not in nc_out.dimensions: + nc_out.__setattr__(a ,nc_in.__getattr__(a)) # Copy variables corresponding to dimensions used by variable copy_vars = list(variable.dimensions) @@ -46,31 +68,50 @@ else: time_bounds_vars = [] for v in copy_vars+time_bounds_vars: if v in nc_in.variables: - hv = nc_out.createVariable(v, nc_in.variables[v].dtype, nc_in.variables[v].dimensions) - for a in nc_in.variables[v].ncattrs(): - hv.setncattr(a, nc_in.variables[v].__getattr__(a)) + if v not in nc_out.variables: + hv = nc_out.createVariable(v, nc_in.variables[v].dtype, nc_in.variables[v].dimensions) + for a in nc_in.variables[v].ncattrs(): + hv.setncattr(a, nc_in.variables[v].__getattr__(a)) + if v in variable.dimensions: + if not nc_in.dimensions[v].isunlimited(): + hv[:] = nc_in.variables[v][:] + print(str(v)) if v in variable.dimensions: - if not nc_in.dimensions[v].isunlimited(): - hv[:] = nc_in.variables[v][:] - else: + if nc_in.dimensions[v].isunlimited(): intime = nc_in.variables[v] - time = hv # Keep around for later + time = intime[:] # Keep around for later # Create new variables -new_mean = nc_out.createVariable(args.variable, variable.dtype, variable.dimensions) -new_squared = nc_out.createVariable(args.variable+'_squared', variable.dtype, variable.dimensions) -new_variance = nc_out.createVariable(args.variable+'_var', variable.dtype, variable.dimensions) +if args.variable not in nc_out.variables: + new_mean = nc_out.createVariable(args.variable, variable.dtype, variable.dimensions) + do_mean = True +else: + do_mean = False + +if args.variable+'_squared' not in nc_out.variables: + new_squared = nc_out.createVariable(args.variable+'_squared', variable.dtype, variable.dimensions) + do_squared = True +else: + do_squared = False + +if args.variable+'_var' not in nc_out.variables: + new_variance = nc_out.createVariable(args.variable+'_var', variable.dtype, variable.dimensions) + do_variance = True +else: + do_variance = False + + for a in variable.ncattrs(): - new_mean.setncattr(a, variable.__getattr__(a)) + if do_mean: new_mean.setncattr(a, variable.__getattr__(a)) if a == 'long_name': - new_squared.setncattr(a, 'Square of '+variable.__getattr__(a)) - new_variance.setncattr(a, 'Variance of '+variable.__getattr__(a)) + if do_squared: new_squared.setncattr(a, 'Square of '+variable.__getattr__(a)) + if do_variance: new_variance.setncattr(a, 'Variance of '+variable.__getattr__(a)) elif a == 'units': - new_squared.setncattr(a, '('+variable.__getattr__(a)+')^2') - new_variance.setncattr(a, '('+variable.__getattr__(a)+')^2') + if do_squared: new_squared.setncattr(a, '('+variable.__getattr__(a)+')^2') + if do_variance: new_variance.setncattr(a, '('+variable.__getattr__(a)+')^2') else: - new_squared.setncattr(a, variable.__getattr__(a)) - new_variance.setncattr(a, variable.__getattr__(a)) + if do_squared: new_squared.setncattr(a, variable.__getattr__(a)) + if do_variance: new_variance.setncattr(a, variable.__getattr__(a)) numpy.seterr(divide='ignore', invalid='ignore', over='ignore') # To avoid warnings record = -1 @@ -97,14 +138,15 @@ count += 1. mean_val = mean_val/count squared_val = squared_val/count - new_mean[month,:] = mean_val - new_squared[month,:] = squared_val - new_variance[month,:] = squared_val - mean_val**2 - time[month] = time_val/count - if len(time_bounds_vars): - nc_out.variables[time_bounds_vars[0]][month] = b0 - nc_out.variables[time_bounds_vars[1]][month] = b1 - nc_out.variables[time_bounds_vars[2]][month] = count + if do_mean: new_mean[month,:] = mean_val + if do_squared: new_squared[month,:] = squared_val + if do_variance: new_variance[month,:] = squared_val - mean_val**2 + if append is False: + time[month] = time_val/count + if len(time_bounds_vars): + nc_out.variables[time_bounds_vars[0]][month] = b0 + nc_out.variables[time_bounds_vars[1]][month] = b1 + nc_out.variables[time_bounds_vars[2]][month] = count nc_out.close() nc_in.close() From c7c337503ebc98feb0aea508fd2670f40440975c Mon Sep 17 00:00:00 2001 From: John Krasting Date: Wed, 4 Oct 2017 08:07:42 -0400 Subject: [PATCH 13/22] Updated diag table for transport diags - T_ady_2d and T_adx_2d now in ocean_month output stream (needed for CMIP6) - Added umo to ocean_month output stream --- ice_ocean_SIS2/OM4_025/diag_table.MOM6 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/ice_ocean_SIS2/OM4_025/diag_table.MOM6 b/ice_ocean_SIS2/OM4_025/diag_table.MOM6 index 97cadbd276..19be1c24a7 100644 --- a/ice_ocean_SIS2/OM4_025/diag_table.MOM6 +++ b/ice_ocean_SIS2/OM4_025/diag_table.MOM6 @@ -162,7 +162,7 @@ "ocean_model_z","vo", "vo", "ocean_month_z", "all", "mean", "none",2 #"ocean_model", "vo", "vo", "ocean_month", "all", "mean", "none",2 #"ocean_model", "umo", "umo", "ocean_annual", "all", "mean", "none",2 - "ocean_model_z","umo", "umo", "ocean_annual_z", "all", "mean", "none",2 + "ocean_model_z","umo", "umo", "ocean_month_z", "all", "mean", "none",2 #"ocean_model_z","umo", "umo", "ocean_month_z", "all", "mean", "none",2 #"ocean_model", "vmo", "vmo", "ocean_annual", "all", "mean", "none",2 "ocean_model_z","vmo", "vmo", "ocean_annual_z", "all", "mean", "none",2 @@ -185,16 +185,16 @@ "ocean_model_z","vh", "vh", "ocean_annual_z", "all", "mean", "none",2 "ocean_model_z","T_adx", "T_adx", "ocean_annual_z", "all", "mean", "none",2 "ocean_model_z","T_ady", "T_ady", "ocean_annual_z", "all", "mean", "none",2 - "ocean_model", "T_adx_2d", "T_adx_2d", "ocean_annual", "all", "mean", "none",2 - "ocean_model", "T_ady_2d", "T_ady_2d", "ocean_annual", "all", "mean", "none",2 + "ocean_model", "T_adx_2d", "T_adx_2d", "ocean_month", "all", "mean", "none",2 + "ocean_model", "T_ady_2d", "T_ady_2d", "ocean_month", "all", "mean", "none",2 "ocean_model_z","S_adx", "S_adx", "ocean_annual_z", "all", "mean", "none",2 "ocean_model_z","S_ady", "S_ady", "ocean_annual_z", "all", "mean", "none",2 - "ocean_model", "S_adx_2d", "S_adx_2d", "ocean_annual", "all", "mean", "none",2 - "ocean_model", "S_ady_2d", "S_ady_2d", "ocean_annual", "all", "mean", "none",2 - "ocean_model", "ndiff_tracer_trans_x_2d_T","ndiff_tracer_trans_x_2d_T", "ocean_annual", "all", "mean", "none",2 - "ocean_model", "ndiff_tracer_trans_y_2d_T","ndiff_tracer_trans_y_2d_T", "ocean_annual", "all", "mean", "none",2 - "ocean_model", "ndiff_tracer_trans_x_2d_S","ndiff_tracer_trans_x_2d_S", "ocean_annual", "all", "mean", "none",2 - "ocean_model", "ndiff_tracer_trans_y_2d_S","ndiff_tracer_trans_y_2d_S", "ocean_annual", "all", "mean", "none",2 + "ocean_model", "S_adx_2d", "S_adx_2d", "ocean_month", "all", "mean", "none",2 + "ocean_model", "S_ady_2d", "S_ady_2d", "ocean_month", "all", "mean", "none",2 + "ocean_model", "ndiff_tracer_trans_x_2d_T","ndiff_tracer_trans_x_2d_T", "ocean_month", "all", "mean", "none",2 + "ocean_model", "ndiff_tracer_trans_y_2d_T","ndiff_tracer_trans_y_2d_T", "ocean_month", "all", "mean", "none",2 + "ocean_model", "ndiff_tracer_trans_x_2d_S","ndiff_tracer_trans_x_2d_S", "ocean_month", "all", "mean", "none",2 + "ocean_model", "ndiff_tracer_trans_y_2d_S","ndiff_tracer_trans_y_2d_S", "ocean_month", "all", "mean", "none",2 # Density space diagnostics (not necessarily using CMIP names but needed to generate CMIP output in post-processing) "ocean_model_rho2","umo", "umo", "ocean_annual_rho2", "all", "mean", "none",2 From f4c5c9525940ab457904cf7d35b5057157ce0554 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Wed, 4 Oct 2017 14:45:10 -0400 Subject: [PATCH 14/22] Added thkcello/volcello for rho2 output --- ice_ocean_SIS2/OM4_025/diag_table.MOM6 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/ice_ocean_SIS2/OM4_025/diag_table.MOM6 b/ice_ocean_SIS2/OM4_025/diag_table.MOM6 index 19be1c24a7..f2600e12f6 100644 --- a/ice_ocean_SIS2/OM4_025/diag_table.MOM6 +++ b/ice_ocean_SIS2/OM4_025/diag_table.MOM6 @@ -71,6 +71,8 @@ # Generally save annuals, and sometimes monthly and daily. "ocean_model_z", "volcello", "volcello", "ocean_annual_z", "all", "mean", "none",2 # Cell measure for 3d data "ocean_model_z", "volcello", "volcello", "ocean_month_z", "all", "mean", "none",2 # Cell measure for 3d data + "ocean_model_rho2", "volcello", "volcello", "ocean_annual_rho2", "all", "mean", "none",2 # Cell measure for 3d data + "ocean_model_rho2", "volcello", "volcello", "ocean_month_rho2", "all", "mean", "none",2 # Cell measure for 3d data "ocean_model", "volcello", "volcello", "ocean_annual", "all", "mean", "none",2 # Cell measure for 3d data #"ocean_model", "volcello", "volcello", "ocean_month", "all", "mean", "none",2 # Cell measure for 3d data "ocean_model", "pbo", "pbo", "ocean_annual", "all", "mean", "none",2 @@ -82,6 +84,8 @@ "ocean_model", "masso", "masso", "ocean_scalar_month", "all", "mean", "none",2 # global sum masscello "ocean_model", "masso", "masso", "ocean_scalar_annual", "all", "mean", "none",2 # global sum masscello "ocean_model", "thkcello", "thkcello", "ocean_annual", "all", "mean", "none",2 # Only needed in native space a static field needs to be provided for CMIP6 + "ocean_model_rho2", "thkcello", "thkcello", "ocean_annual_rho2", "all", "mean", "none",2 # Only needed in native space a static field needs to be provided for CMIP6 + "ocean_model_rho2", "thkcello", "thkcello", "ocean_month_rho2", "all", "mean", "none",2 # Only needed in native space a static field needs to be provided for CMIP6 #"ocean_model", "thkcello", "thkcello", "ocean_month", "all", "mean", "none",2 "ocean_model", "volo", "volo", "ocean_scalar_month", "all", "mean", "none",2 # global sum thkcello "ocean_model", "volo", "volo", "ocean_scalar_annual", "all", "mean", "none",2 # global sum thkcello From ab64e2aa8828a95b8517e3129b9ebe0329931c1a Mon Sep 17 00:00:00 2001 From: John Krasting Date: Wed, 4 Oct 2017 15:16:24 -0400 Subject: [PATCH 15/22] Various fixes for mfo (straits) diagnostic - Added 'runpath' to filenames in strait_transport_timeseries.py - Added additional strlen2 dimension - Added netCDF attributes for mfo and wmo diagnostics - Cast mfo as np.array before posting (netCDF4 Python did not like 0-D mask attributes) --- tools/analysis/refineDiag_ocean_month_z.py | 49 +++++++++++-------- tools/analysis/strait_transport_timeseries.py | 12 ++--- 2 files changed, 35 insertions(+), 26 deletions(-) diff --git a/tools/analysis/refineDiag_ocean_month_z.py b/tools/analysis/refineDiag_ocean_month_z.py index 4d2c0b279a..e4b5496aba 100755 --- a/tools/analysis/refineDiag_ocean_month_z.py +++ b/tools/analysis/refineDiag_ocean_month_z.py @@ -156,15 +156,15 @@ def main(args): wmo.standard_name = 'upward_ocean_mass_transport' #-- mfo - #_, mfo, straits = sum_transport_in_straits(args.straitdir, monthly_average = True) - #mfo[mfo].mask = nc_misval - #mfo = np.ma.array(mfo, fill_value=nc_misval) - #strait_names = np.array( [strait.cmor_name for strait in straits] ) - #mfo.long_name = 'Sea Water Transport' - #mfo.units = 'kg s-1' - #mfo.cell_methods = 'time:mean' - #mfo.time_avg_info = 'average_T1,average_T2,average_DT' - #mfo.standard_name = 'sea_water_transport_across_line' + _, mfo, straits = sum_transport_in_straits(args.straitdir, monthly_average = True) + #mfo[mfo.mask] = nc_misval + mfo = np.ma.array(mfo, fill_value=nc_misval) + strait_names = np.array( [strait.cmor_name for strait in straits] ) + mfo.long_name = 'Sea Water Transport' + mfo.units = 'kg s-1' + mfo.cell_methods = 'time:mean' + mfo.time_avg_info = 'average_T1,average_T2,average_DT' + mfo.standard_name = 'sea_water_transport_across_line' #-- Read time bounds nv = f_in.variables['nv'] @@ -202,8 +202,9 @@ def main(args): time_dim = f_out.createDimension('time', size=None) basin_dim = f_out.createDimension('basin', size=3) - #strait_dim = f_out.createDimension('strait', size=len(straits)) - strlen_dim = f_out.createDimension('strlen', size=21) + strait_dim = f_out.createDimension('strait', size=len(straits)) + strlen1_dim = f_out.createDimension('strlen1', size=21) + strlen2_dim = f_out.createDimension('strlen2', size=31) xh_dim = f_out.createDimension('xh', size=len(xh[:])) yh_dim = f_out.createDimension('yh', size=len(yh[:])) yq_dim = f_out.createDimension('yq', size=len(yq[:])) @@ -212,11 +213,11 @@ def main(args): nv_dim = f_out.createDimension('nv', size=len(nv[:])) time_out = f_out.createVariable('time', np.float64, ('time')) - #basin_out = f_out.createVariable('basin', np.int32, ('basin')) + xh_out = f_out.createVariable('xh', np.float64, ('xh')) yh_out = f_out.createVariable('yh', np.float64, ('yh')) yq_out = f_out.createVariable('yq', np.float64, ('yq')) - region_out = f_out.createVariable('region', 'c', ('basin', 'strlen')) - #strait_out = f_out.createVariable('strait', 'c', ('strait', 'strlen')) + region_out = f_out.createVariable('region', 'c', ('basin', 'strlen1')) + strait_out = f_out.createVariable('strait', 'c', ('strait', 'strlen2')) z_l_out = f_out.createVariable('z_l', np.float64, ('z_l')) z_i_out = f_out.createVariable('z_i', np.float64, ('z_i')) nv_out = f_out.createVariable('nv', np.float64, ('nv')) @@ -233,8 +234,8 @@ def main(args): wmo_out = f_out.createVariable('wmo', np.float32, ('time', 'z_i', 'yh', 'xh'), fill_value=nc_misval) wmo_out.missing_value = nc_misval - #mfo_out = f_out.createVariable('wmo', np.float32, ('time', 'strait'), fill_value=nc_misval) - #mfo_out.missing_value = nc_misval + mfo_out = f_out.createVariable('mfo', np.float32, ('time', 'strait'), fill_value=nc_misval) + mfo_out.missing_value = nc_misval average_T1_out = f_out.createVariable('average_T1', np.float64, ('time')) average_T2_out = f_out.createVariable('average_T2', np.float64, ('time')) @@ -242,7 +243,7 @@ def main(args): time_bnds_out = f_out.createVariable('time_bnds', np.float64, ('time', 'nv')) time_out.setncatts(tax.__dict__) - #xh_out.setncatts(xh.__dict__) + xh_out.setncatts(xh.__dict__) yh_out.setncatts(yh.__dict__) yq_out.setncatts(yq.__dict__) z_l_out.setncatts(z_l.__dict__) @@ -258,13 +259,19 @@ def main(args): for k in msftyzmpa.__dict__.keys(): if k[0] != '_': msftyzmpa_out.setncattr(k,msftyzmpa.__dict__[k]) + for k in mfo.__dict__.keys(): + if k[0] != '_': mfo_out.setncattr(k,mfo.__dict__[k]) + + for k in wmo.__dict__.keys(): + if k[0] != '_': wmo_out.setncattr(k,wmo.__dict__[k]) + average_T1_out.setncatts(average_T1.__dict__) average_T2_out.setncatts(average_T2.__dict__) average_DT_out.setncatts(average_DT.__dict__) time_bnds_out.setncatts(time_bnds.__dict__) time_out[:] = np.array(tax[:]) - #xh_out[:] = np.array(xh[:]) + xh_out[:] = np.array(xh[:]) yh_out[:] = np.array(yh[:]) yq_out[:] = np.array(yq[:]) z_l_out[:] = np.array(z_l[:]) @@ -275,15 +282,17 @@ def main(args): msftyzsmpa_out[:] = np.ma.array(msftyzsmpa[:]) msftyzmpa_out[:] = np.ma.array(msftyzmpa[:]) wmo_out[:] = np.ma.array(wmo[:]) - #mfo_out[:] = np.ma.array(mfo[:]) + mfo_out[:] = np.array(mfo[:]) average_T1_out[:] = average_T1[:] average_T2_out[:] = average_T2[:] average_DT_out[:] = average_DT[:] time_bnds_out[:] = time_bnds[:] + print(region.shape) + print(strait_names.shape) region_out[:] = nc.stringtochar(region) - #strait_out[:] = nc.stringtochar(strait_names) + strait_out[:] = nc.stringtochar(strait_names) f_out.close() exit(0) diff --git a/tools/analysis/strait_transport_timeseries.py b/tools/analysis/strait_transport_timeseries.py index 9117164f42..1fe65698c9 100644 --- a/tools/analysis/strait_transport_timeseries.py +++ b/tools/analysis/strait_transport_timeseries.py @@ -19,20 +19,20 @@ def sum_transport_in_straits(runpath, monthly_average = False): if (len(u_file)==0 and len(v_file)==0): print("Warning: File not found for %s" % strait[sidx].mom6_name) continue - u_vargroup = Dataset(u_file[0]).variables - v_vargroup = Dataset(v_file[0]).variables + u_vargroup = Dataset(runpath+'/'+u_file[0]).variables + v_vargroup = Dataset(runpath+'/'+v_file[0]).variables elif strait[sidx].is_zonal: v_file = fn_filter(listdir(runpath), '*' + strait[sidx].mom6_name + '*.nc') if (len(v_file)==0): print("Warning: File not found for %s" % strait[sidx].mom6_name) continue - v_vargroup = Dataset(v_file[0]).variables + v_vargroup = Dataset(runpath+'/'+v_file[0]).variables elif strait[sidx].is_meridional: u_file = fn_filter(listdir(runpath), '*' + strait[sidx].mom6_name + '*.nc') if (len(u_file)==0): print("Warning: File not found for %s" % strait[sidx].mom6_name) continue - u_vargroup = Dataset(u_file[0]).variables + u_vargroup = Dataset(runpath+'/'+u_file[0]).variables # Need to find the first interface deeper than or equal to the requested z-limit. If deeper, then we'll need to scale the # bottommost part of the column @@ -47,7 +47,7 @@ def sum_transport_in_straits(runpath, monthly_average = False): else: vmo = v_vargroup['vmo'][:,:,:,:] strait[sidx].transport += vmo.sum(axis=(1,2,3)) - Dataset(v_file[0]).close() + Dataset(runpath+'/'+v_file[0]).close() if strait[sidx].is_meridional: strait[sidx].time = u_vargroup['time'][:] @@ -60,7 +60,7 @@ def sum_transport_in_straits(runpath, monthly_average = False): else: umo = u_vargroup['umo'][:,:,:,:] strait[sidx].transport += umo.sum(axis=(1,2,3)) - Dataset(u_file[0]).close() + Dataset(runpath+'/'+u_file[0]).close() if monthly_average: strait[sidx].transport = make_monthly_averages(strait[sidx].transport) strait[sidx].time = make_monthly_averages(strait[sidx].time) From 29c2b13c52c174d6301fb567244f33d3aa5e8b02 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Wed, 4 Oct 2017 19:06:31 -0400 Subject: [PATCH 16/22] Flexible mask tiling based on time dimension --- tools/analysis/refineDiag_ocean_month.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tools/analysis/refineDiag_ocean_month.py b/tools/analysis/refineDiag_ocean_month.py index 4c5db85bb4..9824a0bb06 100755 --- a/tools/analysis/refineDiag_ocean_month.py +++ b/tools/analysis/refineDiag_ocean_month.py @@ -51,7 +51,8 @@ def heat_trans_by_basin(x,mask=None,lat=None,minlat=None): result.mask = varmask if (minlat is not None) and (lat is not None): if len(lat.shape) == 1: - lat = np.expand_dims(lat,0) + lat = np.tile(lat,0) + lat = np.tile(lat,(len(x),1)) result = np.ma.masked_where(np.less(lat,minlat),result) return result From 8fad2c95d149dff7a1bad02d84f273827b92c73e Mon Sep 17 00:00:00 2001 From: John Krasting Date: Wed, 4 Oct 2017 19:07:22 -0400 Subject: [PATCH 17/22] Miscellaneous clean-ups from testing - removed print statements used in debugging - added '-s' option to MOM6_refineDiag.csh - collpased back to single strlen dimension in refineDiag_ocean_month_z.py --- tools/analysis/MOM6_refineDiag.csh | 2 +- tools/analysis/calc_variance.py | 1 - tools/analysis/refineDiag_ocean_month_z.py | 11 ++++------- 3 files changed, 5 insertions(+), 9 deletions(-) diff --git a/tools/analysis/MOM6_refineDiag.csh b/tools/analysis/MOM6_refineDiag.csh index d1fe08d6c4..2a712825b1 100755 --- a/tools/analysis/MOM6_refineDiag.csh +++ b/tools/analysis/MOM6_refineDiag.csh @@ -104,7 +104,7 @@ $script_dir/EddyKineticEnergy.py -g $ocean_static_file -o $out_dir/refineDiag_o echo '==== Offline Diagnostics ====' $script_dir/refineDiag_ocean_month.py -b $basin_codes_file -r refineDiagDir $yr1.ocean_month.nc -$script_dir/refineDiag_ocean_month_z.py -b $basin_codes_file -r refineDiagDir -s $straitdir $yr1.ocean_month_z.nc +$script_dir/refineDiag_ocean_month_z.py -b $basin_codes_file -r refineDiagDir -s ./ $yr1.ocean_month_z.nc $script_dir/refineDiag_ocean_month_rho2.py -b $basin_codes_file -r refineDiagDir $yr1.ocean_month_rho2.nc $script_dir/calc_variance.py zos $yr1.ocean_daily.nc $refineDiagDir/$yr1.ocean_month_refined.nc diff --git a/tools/analysis/calc_variance.py b/tools/analysis/calc_variance.py index 378fb6f890..5f8084c8d4 100755 --- a/tools/analysis/calc_variance.py +++ b/tools/analysis/calc_variance.py @@ -75,7 +75,6 @@ if v in variable.dimensions: if not nc_in.dimensions[v].isunlimited(): hv[:] = nc_in.variables[v][:] - print(str(v)) if v in variable.dimensions: if nc_in.dimensions[v].isunlimited(): intime = nc_in.variables[v] diff --git a/tools/analysis/refineDiag_ocean_month_z.py b/tools/analysis/refineDiag_ocean_month_z.py index e4b5496aba..ab481ff9c0 100755 --- a/tools/analysis/refineDiag_ocean_month_z.py +++ b/tools/analysis/refineDiag_ocean_month_z.py @@ -68,7 +68,7 @@ def main(args): # Note: The Atlantic should include other smaller bays/seas that are # included in the definition used in meridional_overturning.py - region = np.array(['atlantic_arctic_ocean','indian_pacific_ocean','global_ocean']) + region = np.array(['atlantic_arctic_ocean ','indian_pacific_ocean','global_ocean']) #-- Read basin masks f_basin = nc.Dataset(args.basinfile) @@ -203,8 +203,7 @@ def main(args): time_dim = f_out.createDimension('time', size=None) basin_dim = f_out.createDimension('basin', size=3) strait_dim = f_out.createDimension('strait', size=len(straits)) - strlen1_dim = f_out.createDimension('strlen1', size=21) - strlen2_dim = f_out.createDimension('strlen2', size=31) + strlen_dim = f_out.createDimension('strlen', size=31) xh_dim = f_out.createDimension('xh', size=len(xh[:])) yh_dim = f_out.createDimension('yh', size=len(yh[:])) yq_dim = f_out.createDimension('yq', size=len(yq[:])) @@ -216,8 +215,8 @@ def main(args): xh_out = f_out.createVariable('xh', np.float64, ('xh')) yh_out = f_out.createVariable('yh', np.float64, ('yh')) yq_out = f_out.createVariable('yq', np.float64, ('yq')) - region_out = f_out.createVariable('region', 'c', ('basin', 'strlen1')) - strait_out = f_out.createVariable('strait', 'c', ('strait', 'strlen2')) + region_out = f_out.createVariable('region', 'c', ('basin', 'strlen')) + strait_out = f_out.createVariable('strait', 'c', ('strait', 'strlen')) z_l_out = f_out.createVariable('z_l', np.float64, ('z_l')) z_i_out = f_out.createVariable('z_i', np.float64, ('z_i')) nv_out = f_out.createVariable('nv', np.float64, ('nv')) @@ -289,8 +288,6 @@ def main(args): average_DT_out[:] = average_DT[:] time_bnds_out[:] = time_bnds[:] - print(region.shape) - print(strait_names.shape) region_out[:] = nc.stringtochar(region) strait_out[:] = nc.stringtochar(strait_names) f_out.close() From b1c65de249675e6f261154d2d0ee0516781dd669 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Thu, 5 Oct 2017 12:09:31 -0400 Subject: [PATCH 18/22] Updated umo in diag_table - Now in both the monthly and annual streams --- ice_ocean_SIS2/OM4_025/diag_table.MOM6 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ice_ocean_SIS2/OM4_025/diag_table.MOM6 b/ice_ocean_SIS2/OM4_025/diag_table.MOM6 index f2600e12f6..bc14e146c2 100644 --- a/ice_ocean_SIS2/OM4_025/diag_table.MOM6 +++ b/ice_ocean_SIS2/OM4_025/diag_table.MOM6 @@ -166,8 +166,8 @@ "ocean_model_z","vo", "vo", "ocean_month_z", "all", "mean", "none",2 #"ocean_model", "vo", "vo", "ocean_month", "all", "mean", "none",2 #"ocean_model", "umo", "umo", "ocean_annual", "all", "mean", "none",2 - "ocean_model_z","umo", "umo", "ocean_month_z", "all", "mean", "none",2 -#"ocean_model_z","umo", "umo", "ocean_month_z", "all", "mean", "none",2 + "ocean_model_z","umo", "umo", "ocean_annual_z", "all", "mean", "none",2 + "ocean_model_z","umo", "umo", "ocean_month_z", "all", "mean", "none",2 #"ocean_model", "vmo", "vmo", "ocean_annual", "all", "mean", "none",2 "ocean_model_z","vmo", "vmo", "ocean_annual_z", "all", "mean", "none",2 "ocean_model_z","vmo", "vmo", "ocean_month_z", "all", "mean", "none",2 From 9d4536ae48bd14cb5ca3cc9ebc0e8ef8d7470b12 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Thu, 5 Oct 2017 13:10:11 -0400 Subject: [PATCH 19/22] Updated cell_methods attrs. for refined diags --- tools/analysis/refineDiag_ocean_month.py | 2 +- tools/analysis/refineDiag_ocean_month_rho2.py | 4 ++-- tools/analysis/refineDiag_ocean_month_z.py | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/tools/analysis/refineDiag_ocean_month.py b/tools/analysis/refineDiag_ocean_month.py index 9824a0bb06..257bbf3a8d 100755 --- a/tools/analysis/refineDiag_ocean_month.py +++ b/tools/analysis/refineDiag_ocean_month.py @@ -119,7 +119,7 @@ def main(args): hfbasin = np.ma.array(hfbasin,fill_value=nc_misval) hfbasin.long_name = 'Northward Ocean Heat Transport' hfbasin.units = 'W' - hfbasin.cell_methods = 'xh:sum yq:sum basin:mean time:mean' + hfbasin.cell_methods = 'yq:sum time:mean' hfbasin.comment = 'Indo-Pacific heat transport begins at 34 S' hfbasin.time_avg_info = 'average_T1,average_T2,average_DT' hfbasin.standard_name = 'northward_ocean_heat_transport' diff --git a/tools/analysis/refineDiag_ocean_month_rho2.py b/tools/analysis/refineDiag_ocean_month_rho2.py index f9a0122483..23800f70cf 100755 --- a/tools/analysis/refineDiag_ocean_month_rho2.py +++ b/tools/analysis/refineDiag_ocean_month_rho2.py @@ -64,7 +64,7 @@ def main(args): msftyrho = np.ma.array(msftyrho,fill_value=1.e20) msftyrho.long_name = 'Ocean Y Overturning Mass Streamfunction' msftyrho.units = 'kg s-1' - msftyrho.cell_methods = 'rho2_i:sum yq:sum basin:mean time:mean' + msftyrho.cell_methods = 'rho2_i:sum yq:sum time:mean' msftyrho.time_avg_info = 'average_T1,average_T2,average_DT' msftyrho.standard_name = 'ocean_y_overturning_mass_streamfunction' @@ -78,7 +78,7 @@ def main(args): msftyrhompa = np.ma.array(msftyrhompa,fill_value=1.e20) msftyrhompa.long_name = 'ocean Y overturning mass streamfunction due to parameterized mesoscale advection' msftyrhompa.units = 'kg s-1' - msftyrhompa.cell_methods = 'rho2_i:sum yq:sum basin:mean time:mean' + msftyrhompa.cell_methods = 'rho2_i:sum yq:sum time:mean' msftyrhompa.time_avg_info = 'average_T1,average_T2,average_DT' msftyrhompa.standard_name = 'ocean_y_overturning_mass_streamfunction_due_to_parameterized_'+\ 'mesoscale_advection' diff --git a/tools/analysis/refineDiag_ocean_month_z.py b/tools/analysis/refineDiag_ocean_month_z.py index ab481ff9c0..217fdc5583 100755 --- a/tools/analysis/refineDiag_ocean_month_z.py +++ b/tools/analysis/refineDiag_ocean_month_z.py @@ -110,7 +110,7 @@ def main(args): msftyyz = np.ma.array(msftyyz,fill_value=nc_misval) msftyyz.long_name = 'Ocean Y Overturning Mass Streamfunction' msftyyz.units = 'kg s-1' - msftyyz.cell_methods = 'z_i:sum yq:sum basin:mean time:mean' + msftyyz.cell_methods = 'z_i:sum yq:sum time:mean' msftyyz.time_avg_info = 'average_T1,average_T2,average_DT' msftyyz.standard_name = 'ocean_y_overturning_mass_streamfunction' @@ -124,7 +124,7 @@ def main(args): msftyzmpa = np.ma.array(msftyzmpa,fill_value=nc_misval) msftyzmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized mesoscale advection' msftyzmpa.units = 'kg s-1' - msftyzmpa.cell_methods = 'z_i:sum yq:sum basin:mean time:mean' + msftyzmpa.cell_methods = 'z_i:sum yq:sum time:mean' msftyzmpa.time_avg_info = 'average_T1,average_T2,average_DT' msftyzmpa.standard_name = 'ocean_y_overturning_mass_streamfunction_due_to_parameterized_'+\ 'mesoscale_advection' @@ -139,7 +139,7 @@ def main(args): msftyzsmpa = np.ma.array(msftyzsmpa,fill_value=nc_misval) msftyzsmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized submesoscale advection' msftyzsmpa.units = 'kg s-1' - msftyzsmpa.cell_methods = 'z_i:sum yq:sum basin:mean time:mean' + msftyzsmpa.cell_methods = 'z_i:sum yq:sum time:mean' msftyzsmpa.time_avg_info = 'average_T1,average_T2,average_DT' msftyzsmpa.standard_name = 'ocean_meridional_overturning_mass_streamfunction_due_to_parameterized_'+\ 'submesoscale_advection' From 1ebf336b6ea55d43ebd6710aee34d37150057f90 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Thu, 5 Oct 2017 13:11:07 -0400 Subject: [PATCH 20/22] Removed call to calc_varaince.py - Removed call to calc_variance.py from MOM6_refineDiag.csh - 'zos' and 'zos_squared' ('zossq') are now in the diag_table. --- tools/analysis/MOM6_refineDiag.csh | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/tools/analysis/MOM6_refineDiag.csh b/tools/analysis/MOM6_refineDiag.csh index 2a712825b1..5b2fa068b5 100755 --- a/tools/analysis/MOM6_refineDiag.csh +++ b/tools/analysis/MOM6_refineDiag.csh @@ -95,7 +95,6 @@ set ocean_static_file = $yr1.ocean_static.nc if ( -e $yr1.ocean_static_no_mask_table.nc ) set ocean_static_file = $yr1.ocean_static_no_mask_table.nc set basin_codes_file = $yr1.basin_codes.nc -set strait_dir = "./" echo '====annual mean Eddy Kinetic Energy======' mkdir -p $out_dir/refineDiag_ocean_annual/EddyKineticEnergy @@ -103,13 +102,15 @@ $script_dir/EddyKineticEnergy.py -g $ocean_static_file -o $out_dir/refineDiag_o echo '==== Offline Diagnostics ====' -$script_dir/refineDiag_ocean_month.py -b $basin_codes_file -r refineDiagDir $yr1.ocean_month.nc -$script_dir/refineDiag_ocean_month_z.py -b $basin_codes_file -r refineDiagDir -s ./ $yr1.ocean_month_z.nc -$script_dir/refineDiag_ocean_month_rho2.py -b $basin_codes_file -r refineDiagDir $yr1.ocean_month_rho2.nc -$script_dir/calc_variance.py zos $yr1.ocean_daily.nc $refineDiagDir/$yr1.ocean_month_refined.nc +$script_dir/refineDiag_ocean_month.py -b $basin_codes_file -r $refineDiagDir $yr1.ocean_month.nc +$script_dir/refineDiag_ocean_month_z.py -b $basin_codes_file -r $refineDiagDir -s ./ $yr1.ocean_month_z.nc +$script_dir/refineDiag_ocean_month_rho2.py -b $basin_codes_file -r $refineDiagDir $yr1.ocean_month_rho2.nc +#$script_dir/calc_variance.py zos $yr1.ocean_daily.nc $refineDiagDir/$yr1.ocean_month_refined.nc -#-- Note: the calc_variance script should be called LAST, since it appends to the ocean_month_refined.nc file -# that is created first by the other scripts. +#-- Note: The calc_variance script pre-dated refineDiag efforts just prior to the start of the GFDL-CM4 DECK runs. +# Based on the diag_table, it looks like the calc_variance script is no longer needed and is now commented out. +# If it is reactivated, it should be called LAST, since it appends to the ocean_month_refined.nc file that is +# created first by the other scripts. echo " ---------- end yearly analysis ---------- " From 96477439980aa4a9d183399d4b94f154d8687522 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Thu, 5 Oct 2017 13:22:44 -0400 Subject: [PATCH 21/22] Added 'coordinates' attribute to basin/strait vars. --- tools/analysis/refineDiag_ocean_month.py | 3 +++ tools/analysis/refineDiag_ocean_month_rho2.py | 4 ++++ tools/analysis/refineDiag_ocean_month_z.py | 7 +++++++ 3 files changed, 14 insertions(+) diff --git a/tools/analysis/refineDiag_ocean_month.py b/tools/analysis/refineDiag_ocean_month.py index 257bbf3a8d..de7eabffa1 100755 --- a/tools/analysis/refineDiag_ocean_month.py +++ b/tools/analysis/refineDiag_ocean_month.py @@ -119,6 +119,7 @@ def main(args): hfbasin = np.ma.array(hfbasin,fill_value=nc_misval) hfbasin.long_name = 'Northward Ocean Heat Transport' hfbasin.units = 'W' + hfbasin.coordinates = 'region' hfbasin.cell_methods = 'yq:sum time:mean' hfbasin.comment = 'Indo-Pacific heat transport begins at 34 S' hfbasin.time_avg_info = 'average_T1,average_T2,average_DT' @@ -199,6 +200,8 @@ def main(args): for k in hfbasin.__dict__.keys(): if k[0] != '_': hfbasin_out.setncattr(k,hfbasin.__dict__[k]) + region_out.setncattr('standard_name','region') + average_T1_out.setncatts(average_T1.__dict__) average_T2_out.setncatts(average_T2.__dict__) average_DT_out.setncatts(average_DT.__dict__) diff --git a/tools/analysis/refineDiag_ocean_month_rho2.py b/tools/analysis/refineDiag_ocean_month_rho2.py index 23800f70cf..ec715b332c 100755 --- a/tools/analysis/refineDiag_ocean_month_rho2.py +++ b/tools/analysis/refineDiag_ocean_month_rho2.py @@ -64,6 +64,7 @@ def main(args): msftyrho = np.ma.array(msftyrho,fill_value=1.e20) msftyrho.long_name = 'Ocean Y Overturning Mass Streamfunction' msftyrho.units = 'kg s-1' + msftyrho.coordinates = 'region' msftyrho.cell_methods = 'rho2_i:sum yq:sum time:mean' msftyrho.time_avg_info = 'average_T1,average_T2,average_DT' msftyrho.standard_name = 'ocean_y_overturning_mass_streamfunction' @@ -78,6 +79,7 @@ def main(args): msftyrhompa = np.ma.array(msftyrhompa,fill_value=1.e20) msftyrhompa.long_name = 'ocean Y overturning mass streamfunction due to parameterized mesoscale advection' msftyrhompa.units = 'kg s-1' + msftyrhompa.coordinates = 'region' msftyrhompa.cell_methods = 'rho2_i:sum yq:sum time:mean' msftyrhompa.time_avg_info = 'average_T1,average_T2,average_DT' msftyrhompa.standard_name = 'ocean_y_overturning_mass_streamfunction_due_to_parameterized_'+\ @@ -138,6 +140,8 @@ def main(args): msftyrhompa_out = f_out.createVariable('msftyrhompa', np.float32, ('time', 'basin', 'rho2_i', 'yq'), fill_value=1.e20) msftyrhompa_out.missing_value = 1.e20 + + region_out.setncattr('standard_name','region') average_T1_out = f_out.createVariable('average_T1', np.float64, ('time')) average_T2_out = f_out.createVariable('average_T2', np.float64, ('time')) diff --git a/tools/analysis/refineDiag_ocean_month_z.py b/tools/analysis/refineDiag_ocean_month_z.py index 217fdc5583..94fc119237 100755 --- a/tools/analysis/refineDiag_ocean_month_z.py +++ b/tools/analysis/refineDiag_ocean_month_z.py @@ -110,6 +110,7 @@ def main(args): msftyyz = np.ma.array(msftyyz,fill_value=nc_misval) msftyyz.long_name = 'Ocean Y Overturning Mass Streamfunction' msftyyz.units = 'kg s-1' + msftyyz.coordinates = 'region' msftyyz.cell_methods = 'z_i:sum yq:sum time:mean' msftyyz.time_avg_info = 'average_T1,average_T2,average_DT' msftyyz.standard_name = 'ocean_y_overturning_mass_streamfunction' @@ -124,6 +125,7 @@ def main(args): msftyzmpa = np.ma.array(msftyzmpa,fill_value=nc_misval) msftyzmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized mesoscale advection' msftyzmpa.units = 'kg s-1' + msftyzmpa.coordinates = 'region' msftyzmpa.cell_methods = 'z_i:sum yq:sum time:mean' msftyzmpa.time_avg_info = 'average_T1,average_T2,average_DT' msftyzmpa.standard_name = 'ocean_y_overturning_mass_streamfunction_due_to_parameterized_'+\ @@ -139,6 +141,7 @@ def main(args): msftyzsmpa = np.ma.array(msftyzsmpa,fill_value=nc_misval) msftyzsmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized submesoscale advection' msftyzsmpa.units = 'kg s-1' + msftyzsmpa.coordinates = 'region' msftyzsmpa.cell_methods = 'z_i:sum yq:sum time:mean' msftyzsmpa.time_avg_info = 'average_T1,average_T2,average_DT' msftyzsmpa.standard_name = 'ocean_meridional_overturning_mass_streamfunction_due_to_parameterized_'+\ @@ -162,6 +165,7 @@ def main(args): strait_names = np.array( [strait.cmor_name for strait in straits] ) mfo.long_name = 'Sea Water Transport' mfo.units = 'kg s-1' + mfo.coordinates = 'strait' mfo.cell_methods = 'time:mean' mfo.time_avg_info = 'average_T1,average_T2,average_DT' mfo.standard_name = 'sea_water_transport_across_line' @@ -264,6 +268,9 @@ def main(args): for k in wmo.__dict__.keys(): if k[0] != '_': wmo_out.setncattr(k,wmo.__dict__[k]) + region_out.setncattr('standard_name','region') + strait_out.setncattr('standard_name','region') + average_T1_out.setncatts(average_T1.__dict__) average_T2_out.setncatts(average_T2.__dict__) average_DT_out.setncatts(average_DT.__dict__) From bcd2b908b3677b63566e01ea623e18ec485bf6e6 Mon Sep 17 00:00:00 2001 From: John Krasting Date: Thu, 5 Oct 2017 14:54:09 -0400 Subject: [PATCH 22/22] Tidy-up odds and ends after QC - corrected cell_methods - cleaned up global attributes --- tools/analysis/refineDiag_ocean_month.py | 7 ++++--- tools/analysis/refineDiag_ocean_month_rho2.py | 9 +++++---- tools/analysis/refineDiag_ocean_month_z.py | 11 ++++++----- 3 files changed, 15 insertions(+), 12 deletions(-) diff --git a/tools/analysis/refineDiag_ocean_month.py b/tools/analysis/refineDiag_ocean_month.py index de7eabffa1..741e2b9924 100755 --- a/tools/analysis/refineDiag_ocean_month.py +++ b/tools/analysis/refineDiag_ocean_month.py @@ -120,7 +120,7 @@ def main(args): hfbasin.long_name = 'Northward Ocean Heat Transport' hfbasin.units = 'W' hfbasin.coordinates = 'region' - hfbasin.cell_methods = 'yq:sum time:mean' + hfbasin.cell_methods = 'yq:point time:mean' hfbasin.comment = 'Indo-Pacific heat transport begins at 34 S' hfbasin.time_avg_info = 'average_T1,average_T2,average_DT' hfbasin.standard_name = 'northward_ocean_heat_transport' @@ -157,8 +157,9 @@ def main(args): f_out = nc.Dataset(args.outfile, 'w', format='NETCDF3_CLASSIC') f_out.setncatts(f_in.__dict__) - f_out.filename = args.outfile - + f_out.filename = os.path.basename(args.outfile) + f_out.delncattr('associated_files') # not needed for these fields + time_dim = f_out.createDimension('time', size=None) basin_dim = f_out.createDimension('basin', size=3) strlen_dim = f_out.createDimension('strlen', size=21) diff --git a/tools/analysis/refineDiag_ocean_month_rho2.py b/tools/analysis/refineDiag_ocean_month_rho2.py index ec715b332c..84002ec12b 100755 --- a/tools/analysis/refineDiag_ocean_month_rho2.py +++ b/tools/analysis/refineDiag_ocean_month_rho2.py @@ -65,7 +65,7 @@ def main(args): msftyrho.long_name = 'Ocean Y Overturning Mass Streamfunction' msftyrho.units = 'kg s-1' msftyrho.coordinates = 'region' - msftyrho.cell_methods = 'rho2_i:sum yq:sum time:mean' + msftyrho.cell_methods = 'rho2_i:point yq:point time:mean' msftyrho.time_avg_info = 'average_T1,average_T2,average_DT' msftyrho.standard_name = 'ocean_y_overturning_mass_streamfunction' @@ -80,7 +80,7 @@ def main(args): msftyrhompa.long_name = 'ocean Y overturning mass streamfunction due to parameterized mesoscale advection' msftyrhompa.units = 'kg s-1' msftyrhompa.coordinates = 'region' - msftyrhompa.cell_methods = 'rho2_i:sum yq:sum time:mean' + msftyrhompa.cell_methods = 'rho2_i:point yq:point time:mean' msftyrhompa.time_avg_info = 'average_T1,average_T2,average_DT' msftyrhompa.standard_name = 'ocean_y_overturning_mass_streamfunction_due_to_parameterized_'+\ 'mesoscale_advection' @@ -118,8 +118,9 @@ def main(args): f_out = nc.Dataset(args.outfile, 'w', format='NETCDF3_CLASSIC') f_out.setncatts(f_in.__dict__) - f_out.filename = args.outfile - + f_out.filename = os.path.basename(args.outfile) + f_out.delncattr('associated_files') # not needed for these fields + time_dim = f_out.createDimension('time', size=None) basin_dim = f_out.createDimension('basin', size=3) strlen_dim = f_out.createDimension('strlen', size=21) diff --git a/tools/analysis/refineDiag_ocean_month_z.py b/tools/analysis/refineDiag_ocean_month_z.py index 94fc119237..609c7a9b4e 100755 --- a/tools/analysis/refineDiag_ocean_month_z.py +++ b/tools/analysis/refineDiag_ocean_month_z.py @@ -111,7 +111,7 @@ def main(args): msftyyz.long_name = 'Ocean Y Overturning Mass Streamfunction' msftyyz.units = 'kg s-1' msftyyz.coordinates = 'region' - msftyyz.cell_methods = 'z_i:sum yq:sum time:mean' + msftyyz.cell_methods = 'z_i:point yq:point time:mean' msftyyz.time_avg_info = 'average_T1,average_T2,average_DT' msftyyz.standard_name = 'ocean_y_overturning_mass_streamfunction' @@ -126,7 +126,7 @@ def main(args): msftyzmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized mesoscale advection' msftyzmpa.units = 'kg s-1' msftyzmpa.coordinates = 'region' - msftyzmpa.cell_methods = 'z_i:sum yq:sum time:mean' + msftyzmpa.cell_methods = 'z_i:point yq:point time:mean' msftyzmpa.time_avg_info = 'average_T1,average_T2,average_DT' msftyzmpa.standard_name = 'ocean_y_overturning_mass_streamfunction_due_to_parameterized_'+\ 'mesoscale_advection' @@ -142,7 +142,7 @@ def main(args): msftyzsmpa.long_name = 'ocean Y overturning mass streamfunction due to parameterized submesoscale advection' msftyzsmpa.units = 'kg s-1' msftyzsmpa.coordinates = 'region' - msftyzsmpa.cell_methods = 'z_i:sum yq:sum time:mean' + msftyzsmpa.cell_methods = 'z_i:point yq:point time:mean' msftyzsmpa.time_avg_info = 'average_T1,average_T2,average_DT' msftyzsmpa.standard_name = 'ocean_meridional_overturning_mass_streamfunction_due_to_parameterized_'+\ 'submesoscale_advection' @@ -154,9 +154,10 @@ def main(args): wmo = np.ma.array(wmo,fill_value=nc_misval) wmo.long_name = 'Upward mass transport from resolved and parameterized advective transport' wmo.units = 'kg s-1' - wmo.cell_methods = 'z_i:sum xh:sum yh:sum time:mean' + wmo.cell_methods = 'z_i:point xh:sum yh:sum time:mean' wmo.time_avg_info = 'average_T1,average_T2,average_DT' wmo.standard_name = 'upward_ocean_mass_transport' + wmo.cell_measures = 'area:areacello' #-- mfo _, mfo, straits = sum_transport_in_straits(args.straitdir, monthly_average = True) @@ -202,7 +203,7 @@ def main(args): f_out = nc.Dataset(args.outfile, 'w', format='NETCDF3_CLASSIC') f_out.setncatts(f_in.__dict__) - f_out.filename = args.outfile + f_out.filename = os.path.basename(args.outfile) time_dim = f_out.createDimension('time', size=None) basin_dim = f_out.createDimension('basin', size=3)