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..c94cd0b873 --- /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,:,:] += 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 + 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() +