Skip to content

Commit

Permalink
Merge pull request mom-ocean#163 from jkrasting/basin_refinediag
Browse files Browse the repository at this point in the history 
CMIP6 Refine Diag Script

- Closes mom-ocean#178 .
  • Loading branch information
@adcroft
adcroft committed Oct 5, 2017
2 parents 8d4de1e + bcd2b90 commit f15b916
Show file tree
Hide file tree
Showing 9 changed files with 1,079 additions and 58 deletions.
22 changes: 13 additions & 9 deletions ice_ocean_SIS2/OM4_025/diag_table.MOM6
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -163,7 +167,7 @@
#"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
"ocean_model_z","vmo", "vmo", "ocean_month_z", "all", "mean", "none",2
Expand All @@ -185,16 +189,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
Expand Down
50 changes: 31 additions & 19 deletions tools/analysis/MOM6_refineDiag.csh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <refineDiag> tag of the xmls, e.g.,
# To make this happen a path to the (would be) script should appear
# in the <refineDiag> tag of the xmls, e.g.,
# <refineDiag script="$(NB_ROOT)/mom6/tools/analysis/MOM6_refineDiag.csh"/>
# Note that the above script should exist when frepp is called.
# This could be achieved by cloning the mom6 git repo in the <csh> 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 <csh> section of the setup block
# in the corresponding the gfdl platfrom. E.g.,
# <csh><![CDATA[
# source $MODULESHOME/init/csh
# module use -a /home/John.Krasting/local/modulefiles
# module purge
# module load jpk-analysis/0.0.4
# module load $(FRE_VERSION)
### The following clones the mom6 git repo which should contain all the pp scripts
### The following clones the mom6 git repo which should contain all the pp scripts
# setenv NBROOT /nbhome/$USER/$(FRE_STEM)$(DEBUGLEVEL)
# mkdir -p $NBROOT
# cd $NBROOT
Expand All @@ -45,19 +45,19 @@ echo ""
#Niki: Note that here we do not have any FRE environment to /nbhome
# NBROOT should be set as an environ vatiable at the setup in the xml
# setenv NBROOT /nbhome/$USER/$(FRE_STEM)$(DEBUGLEVEL)
set src_dir=$NBROOT/mom6/tools/analysis
set src_dir=$NBROOT/mom6/tools/analysis
# The following variables are set by frepp, but frepp is not called yet at refineDiag stage of FRE workflow,
# so we need to explicitly set them here
set descriptor = $name
set descriptor = $name
set out_dir = $NBROOT #Niki: How can we set this to frepp analysisdir /nbhome
set yr1 = $oname
set yr2 = $oname
set databegyr = $oname
set dataendyr = $oname
set yr1 = $oname
set yr2 = $oname
set databegyr = $oname
set dataendyr = $oname
set datachunk = 1
#Try setting fre version to the caller version
if ( ! $?FREVERSION ) set FREVERSION = fre
set fremodule = $FREVERSION
set fremodule = $FREVERSION
set freanalysismodule = fre-analysis/test

# make sure valid platform and required modules are loaded
Expand All @@ -82,7 +82,7 @@ if (! $?FRE_ANALYSIS_HOME) then
endif

#
#At this point of the FRE workflow we are in a /vftmp directory on an analysis node
#At this point of the FRE workflow we are in a /vftmp directory on an analysis node
#with all the history files for the current finished year ${yr1} unpacked and present
#
echo "We are inside the refineDiag script"
Expand All @@ -91,14 +91,26 @@ 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


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

#-- 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 ---------- "

Expand Down
101 changes: 71 additions & 30 deletions tools/analysis/calc_variance.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.')
Expand All @@ -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]
Expand All @@ -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)
Expand All @@ -46,31 +68,49 @@
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][:]
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
Expand All @@ -97,14 +137,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()
Expand Down
30 changes: 30 additions & 0 deletions tools/analysis/m6toolbox.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,6 +169,36 @@ 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 = 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) * -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):
"""
Find (j,i) of cell with center nearest to (x0,y0).
Expand Down
Loading

0 comments on commit f15b916

Please sign in to comment.