CMIP6 Refine Diag Script

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
@@ -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
@@ -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

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

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,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
@@ -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()

tools/analysis/m6toolbox.py

@@ -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).