Add global moc postprocess

config.analysis

@@ -82,6 +82,9 @@ generate = 1
 [seaice_modelvsobs]
 generate = 1
 
+[moc_postprocess]
+generate = 0
+
 [ohc_timeseries]
 ## compare to output from another model run?
 #compare_with_model = True
@@ -131,6 +134,18 @@ compare_with_obs = True
 # output, N_movavg=12 corresponds to a 12-month moving average window)
 N_movavg = 12
 
+
+[moc_postprocess]
+filename  = ../mpaso.hist.am.timeSeriesStats.yr100-119.avg.nc
+filenamemesh = ../mpaso.rst.0002-01-01_00000.nc
+filenamemask = ../EC60to30v3_SingleRegionAtlanticWTransportTransects_masks.nc
+timeIndex = 0
+horVelName = time_avg_normalVelocity
+vertVelName = time_avg_vertVelocityTop
+contourLines = np.arange(-8,20.1,2)
+# This is specific to Atlantic MOC:
+mocLat = np.arange(-34.5,70.1,0.5)
+
 [seaice_timeseries]
 ## compare to output from another model run?
 #compare_with_model = True

mpas_analysis/configuration/MpasAnalysisConfigParser.py

@@ -1,16 +1,17 @@
-import numbers
-import ast
-
 """
 A configuratin parser class for MPAS analysis.  MpasAnalysisConfigParser adds
 the capabilities to get an option including a default value
 (`getWithDefault(section, option, default, ...)`) and to get options
 that are lists, tuples, dicts, etc (`getExpression(section, option)`).
 
-Author: Xylar Asay-Davis
-Last Modified: 12/07/2016
+Author: Xylar Asay-Davis, Phillip J. Wolfram
+Last Modified: 01/31/2017
 """
 
+import numbers
+import ast
+import numpy as np
+import numpy
 from ConfigParser import ConfigParser
 
 
@@ -44,7 +45,7 @@ def getWithDefault(self, section, option, default):
         self.set(section, option, str(default))
         return default
 
-    def getExpression(self, section, option, elementType=None):
+    def getExpression(self, section, option, elementType=None, usenumpy=False):
         """
         Get an option as an expression (typically a list, though tuples and
         dicts should also work).  `section` and `option` work as in `get(...)`.
@@ -56,11 +57,17 @@ def getExpression(self, section, option, elementType=None):
         useful for ensuring that all elements in a list of numbers are of type
         float, rather than int, when the distinction is important.
 
-        Author: Xylar Asay-Davis
-        Last Modified: 12/0y/2016
+        If `numpy` is True, expression is evaluated within the context of
+        having numpy and / or np available.
+
+        Author: Xylar Asay-Davis, Phillip J. Wolfram
+        Last Modified: 01/31/2017
         """
         expressionString = self.get(section, option)
-        result =  ast.literal_eval(expressionString)
+        if usenumpy:
+            result = eval(expressionString)
+        else:
+            result =  ast.literal_eval(expressionString)
 
         if elementType is not None:
             if isinstance(result, (list, tuple)):

mpas_analysis/ocean/meridional_overturning_circulation.py

@@ -0,0 +1,128 @@
+#!/usr/bin/env python
+"""
+Compute and plot moc, as a post-processing computation
+
+Author: Mark Petersen, Phillip J. Wolfram
+Last Modified: 01/20/2017
+"""
+
+from netCDF4 import Dataset
+import numpy as np
+from ..shared.constants.constants import m3ps_to_Sv
+from ..shared.plot.plotting import plot_moc
+
+def moc_postprocess(config, streamMap=None, variableMap=None): #{{{
+    """
+    Plots a comparison of ACME/MPAS output to SST or MLD observations
+
+    config is an instance of MpasAnalysisConfigParser containing configuration
+    options.
+
+    If present, streamMap is a dictionary of MPAS-O stream names that map to
+    their mpas_analysis counterparts.
+
+    If present, variableMap is a dictionary of MPAS-O variable names that map
+    to their mpas_analysis counterparts.
+
+    Authors: Mark Petersen, Phillip J. Wolfram
+    Modified: 01/20/2017
+    """
+
+    # read parameters from config file
+    fileName = config.get('moc_postprocess', 'filename')
+    fileNameMesh = config.get('moc_postprocess', 'filenamemesh')
+    fileNameMask = config.get('moc_postprocess', 'filenamemask')
+    timeIndex = config.get('moc_postprocess', 'timeIndex')
+    horVelName = config.get('moc_postprocess', 'horVelName')
+    vertVelName = config.get('moc_postprocess', 'vertVelname')
+    mocLat = config.getExpression('moc_postprocess', 'mocLat', usenumpy=True)
+
+    # open files, load variables
+    print('** load variables: \n')
+
+    print 'Reading data from: ' + fileName
+    ncfile = Dataset(fileName,'r')
+    nCells = ncfile.dimensions['nCells'].size
+    nVertLevels = ncfile.dimensions['nVertLevels'].size
+
+    horVel = ncfile.variables[horVelName][:,:,:]
+    vertVel = ncfile.variables[vertVelName][:,:,:]
+
+    print 'Reading data from: ' + fileNameMesh
+    ncfileMesh = Dataset(fileNameMesh,'r')
+    dvEdge = ncfileMesh.variables['dvEdge'][:]
+    areaCell = ncfileMesh.variables['areaCell'][:]
+    latCell = ncfileMesh.variables['latCell'][:]*180./np.pi
+    refBottomDepth = ncfileMesh.variables['refBottomDepth'][:]
+
+    refLayerThickness = np.zeros(nVertLevels)
+    refLayerThickness[0] = refBottomDepth[0]
+    refLayerThickness[1:nVertLevels] = refBottomDepth[1:nVertLevels]-refBottomDepth[0:nVertLevels-1]
+
+    refTopDepth = np.zeros(nVertLevels+1)
+    refTopDepth[1:nVertLevels+1] = refBottomDepth[0:nVertLevels]
+
+    print 'Reading transects from: ' + fileNameMask
+    ncfileMask = Dataset(fileNameMask,'r')
+    nTransects = ncfileMask.dimensions['nTransects'].size
+    maxEdgesInTransect = ncfileMask.dimensions['maxEdgesInTransect'].size
+
+    transectEdgeMaskSigns = ncfileMask.variables['transectEdgeMaskSigns'][:,:]
+    transectEdgeGlobalIDs = ncfileMask.variables['transectEdgeGlobalIDs'][:,:]
+
+    print 'Reading regionss from: ' + fileNameMask
+    ncfileMask = Dataset(fileNameMask,'r')
+    nRegions = ncfileMask.dimensions['nRegions'].size
+    regionCellMasks = ncfileMask.variables['regionCellMasks']
+
+    # compute transport through transects
+    print('** compute transport: \n')
+
+    # the volume transport
+    transport = np.zeros(nTransects);
+    transportZ = np.zeros([nVertLevels,nTransects]);
+    transportZEdge = np.zeros([nVertLevels,maxEdgesInTransect,nTransects]);
+
+    for iTransect in range(nTransects):
+        for i in range(maxEdgesInTransect):
+            if transectEdgeGlobalIDs[iTransect, i]==0:
+                break
+            iEdge = transectEdgeGlobalIDs[iTransect, i] - 1 # subtract 1 because of python 0-indexing
+            for k in range(nVertLevels):
+                transportZEdge[k,i,iTransect] = ( transectEdgeMaskSigns[iEdge,iTransect]
+                        * horVel[timeIndex,iEdge,k] * dvEdge[iEdge] * refLayerThickness[k]*m3ps_to_Sv )
+                transportZ[k,iTransect] = transportZ[k,iTransect] + transportZEdge[k,i,iTransect]
+                transport[iTransect] = transport[iTransect] + transportZEdge[k,i,iTransect]
+
+    # compute MOC
+    print('** compute moc: \n')
+
+    nLat = np.size(mocLat)
+    mocTop = np.zeros((nLat,nVertLevels+1))
+
+    for iRegion in range(nRegions):
+        # assume transects and regions have the same index ordering:
+       iTransect = iRegion
+
+       for k in range(1,nVertLevels+1):
+           mocTop[0,k] = mocTop[0,k-1] + transportZ[k-1,iTransect]/m3ps_to_Sv
+
+       for iLat in range(1,nLat):
+           ind =np.logical_and(np.logical_and(regionCellMasks[:,iRegion]==1,latCell >= mocLat[iLat-1] ), latCell <  mocLat[iLat])
+
+           for k in range(1,nVertLevels+1):
+               mocTop[iLat,k] = mocTop[iLat-1,k] + np.sum(vertVel[timeIndex,ind,k]*areaCell[ind])
+
+       # convert m^3/s to Sverdrup
+       mocTop = mocTop * m3ps_to_Sv
+
+    # plot MOC
+    print('** plot moc: \n')
+    plot_moc(config, mocLat, refTopDepth, mocTop)
+
+    ncfile.close()
+    ncfileMask.close()
+
+    return #}}}
+
+# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python

mpas_analysis/shared/constants/constants.py

@@ -22,4 +22,5 @@
 
 dinmonth = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
 
+m3ps_to_Sv = 1e-6; # m^3/sec flux to Sverdrups