Source code for mpas_analysis.ocean.streamfunction_moc

# -*- coding: utf-8 -*-
# This software is open source software available under the BSD-3 license.
#
# Copyright (c) 2020 Triad National Security, LLC. All rights reserved.
# Copyright (c) 2020 Lawrence Livermore National Security, LLC. All rights
# reserved.
# Copyright (c) 2020 UT-Battelle, LLC. All rights reserved.
#
# Additional copyright and license information can be found in the LICENSE file
# distributed with this code, or at
# https://raw.githubusercontent.com/MPAS-Dev/MPAS-Analysis/master/LICENSE
#

from __future__ import absolute_import, division, print_function, \
    unicode_literals

import xarray as xr
import numpy as np
import netCDF4
import os
from geometric_features import GeometricFeatures
from mpas_tools.ocean.moc import make_moc_basins_and_transects

from mpas_analysis.shared.constants.constants import m3ps_to_Sv
from mpas_analysis.shared.plot import plot_vertical_section_comparison, \
    timeseries_analysis_plot, savefig

from mpas_analysis.shared.io.utility import build_config_full_path, \
    make_directories, get_files_year_month, get_region_mask

from mpas_analysis.shared.io import open_mpas_dataset, write_netcdf

from mpas_analysis.shared.timekeeping.utility import days_to_datetime

from mpas_analysis.shared import AnalysisTask

from mpas_analysis.shared.html import write_image_xml
from mpas_analysis.shared.climatology.climatology import \
    get_climatology_op_directory


[docs]class StreamfunctionMOC(AnalysisTask): # {{{ """ Computation and plotting of model meridional overturning circulation. Will eventually support: * MOC streamfunction, post-processed * MOC streamfunction, from MOC analysis member * MOC time series (max value at 24.5N), post-processed * MOC time series (max value at 24.5N), from MOC analysis member """ # Authors # ------- # Milena Veneziani, Mark Petersen, Phillip Wolfram, Xylar Asay-Davis
[docs] def __init__(self, config, mpasClimatologyTask, controlConfig=None): # {{{ """ Construct the analysis task. Parameters ---------- config : instance of MpasAnalysisConfigParser Contains configuration options mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted controlConfig : ``MpasAnalysisConfigParser``, optional Configuration options for a control run (if any) """ # Authors # ------- # Xylar Asay-Davis # first, call the constructor from the base class (AnalysisTask) super(StreamfunctionMOC, self).__init__( config=config, taskName='streamfunctionMOC', componentName='ocean', tags=['streamfunction', 'moc', 'climatology', 'timeSeries', 'publicObs']) maskSubtask = ComputeMOCMasksSubtask(self) self.add_subtask(maskSubtask) computeClimSubtask = ComputeMOCClimatologySubtask( self, mpasClimatologyTask) computeClimSubtask.run_after(maskSubtask) plotClimSubtask = PlotMOCClimatologySubtask(self, controlConfig) plotClimSubtask.run_after(computeClimSubtask) startYear = config.getint('timeSeries', 'startYear') endYear = config.get('timeSeries', 'endYear') if endYear == 'end': # a valid end year wasn't found, so likely the run was not found, # perhaps because we're just listing analysis tasks endYear = startYear else: endYear = int(endYear) years = range(startYear, endYear + 1) # in the end, we'll combine all the time series into one, but we create # this task first so it's easier to tell it to run after all the # compute tasks combineTimeSeriesSubtask = CombineMOCTimeSeriesSubtask( self, startYears=years, endYears=years) # run one subtask per year for year in years: computeTimeSeriesSubtask = ComputeMOCTimeSeriesSubtask( self, startYear=year, endYear=year) computeTimeSeriesSubtask.run_after(maskSubtask) combineTimeSeriesSubtask.run_after(computeTimeSeriesSubtask) plotTimeSeriesSubtask = PlotMOCTimeSeriesSubtask(self, controlConfig) plotTimeSeriesSubtask.run_after(combineTimeSeriesSubtask)
# }}} # }}} class ComputeMOCMasksSubtask(AnalysisTask): # {{{ """ An analysis subtasks for computing cell masks and southern transects for MOC regions """ # Authors # ------- # Xylar Asay-Davis def __init__(self, parentTask): # {{{ """ Construct the analysis task and adds it as a subtask of the ``parentTask``. Parameters ---------- parentTask : ``AnalysisTask`` The parent task, used to get the ``taskName``, ``config`` and ``componentName`` """ # Authors # ------- # Xylar Asay-Davis # call the constructor from the base class (AnalysisTask) super(ComputeMOCMasksSubtask, self).__init__( config=parentTask.config, taskName=parentTask.taskName, subtaskName='computeRegionMasks', componentName=parentTask.componentName, tags=[]) meshName = self.config.get('input', 'mpasMeshName') self.geojsonFileName = get_region_mask(parentTask.config, 'moc_basins.geojson') self.maskFileName = get_region_mask( self.config, '{}_moc_masks.nc'.format(meshName)) self.maskAndTransectFileName = get_region_mask( self.config, '{}_moc_masks_and_transects.nc'.format(meshName)) # }}} def run_task(self): # {{{ """ Compute the requested climatologies """ # Authors # ------- # Xylar Asay-Davis if os.path.exists(self.maskAndTransectFileName): return config = self.config # make the geojson file gf = GeometricFeatures() mesh_filename = self.runStreams.readpath('restart')[0] maskSubdirectory = build_config_full_path(config, 'output', 'maskSubdirectory') make_directories(maskSubdirectory) make_moc_basins_and_transects(gf, mesh_filename, self.maskAndTransectFileName, geojson_filename=self.geojsonFileName, mask_filename=self.maskFileName, logger=self.logger, dir=maskSubdirectory) # }}} # }}} class ComputeMOCClimatologySubtask(AnalysisTask): # {{{ """ Computation of a climatology of the model meridional overturning circulation. Attributes ---------- mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted """ # Authors # ------- # Milena Veneziani, Mark Petersen, Phillip Wolfram, Xylar Asay-Davis def __init__(self, parentTask, mpasClimatologyTask): # {{{ """ Construct the analysis task. Parameters ---------- parentTask : ``StreamfunctionMOC`` The main task of which this is a subtask mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted """ # Authors # ------- # Xylar Asay-Davis # first, call the constructor from the base class (AnalysisTask) super(ComputeMOCClimatologySubtask, self).__init__( config=parentTask.config, taskName=parentTask.taskName, componentName=parentTask.componentName, tags=parentTask.tags, subtaskName='computeMOCClimatology') self.mpasClimatologyTask = mpasClimatologyTask self.run_after(mpasClimatologyTask) parentTask.add_subtask(self) # }}} def setup_and_check(self): # {{{ """ Perform steps to set up the analysis and check for errors in the setup. Raises ------ ValueError if timeSeriesStatsMonthly is not enabled in the MPAS run """ # Authors # ------- # Xylar Asay-Davis # first, call setup_and_check from the base class (AnalysisTask), # which will perform some common setup, including storing: # self.runDirectory , self.historyDirectory, self.plotsDirectory, # self.namelist, self.runStreams, self.historyStreams, # self.calendar super(ComputeMOCClimatologySubtask, self).setup_and_check() self.startYear = self.mpasClimatologyTask.startYear self.startDate = self.mpasClimatologyTask.startDate self.endYear = self.mpasClimatologyTask.endYear self.endDate = self.mpasClimatologyTask.endDate config = self.config self.mocAnalysisMemberEnabled = self.check_analysis_enabled( analysisOptionName='config_am_mocstreamfunction_enable', raiseException=False) self.sectionName = 'streamfunctionMOC' self.usePostprocessing = config.getExpression( self.sectionName, 'usePostprocessingScript') if not self.usePostprocessing and self.mocAnalysisMemberEnabled: variableList = \ ['timeMonthly_avg_mocStreamvalLatAndDepth', 'timeMonthly_avg_mocStreamvalLatAndDepthRegion'] else: variableList = ['timeMonthly_avg_normalVelocity', 'timeMonthly_avg_vertVelocityTop'] # Add the bolus velocity if GM is enabled try: # the new name self.includeBolus = self.namelist.getbool('config_use_gm') except KeyError: # the old name self.includeBolus = self.namelist.getbool( 'config_use_standardgm') if self.includeBolus: variableList.extend( ['timeMonthly_avg_normalGMBolusVelocity', 'timeMonthly_avg_vertGMBolusVelocityTop']) self.mpasClimatologyTask.add_variables(variableList=variableList, seasons=['ANN']) # }}} def run_task(self): # {{{ """ Process MOC analysis member data if available, or compute MOC at post-processing if not. """ # Authors # ------- # Milena Veneziani, Mark Petersen, Phillip J. Wolfram, Xylar Asay-Davis self.logger.info("Computing climatology of Meridional Overturning " "Circulation (MOC)...") # **** Compute MOC **** if not self.usePostprocessing and self.mocAnalysisMemberEnabled: self._compute_moc_climo_analysismember() else: self._compute_moc_climo_postprocess() # }}} def _compute_moc_climo_analysismember(self): # {{{ """compute mean MOC streamfunction from analysis member""" config = self.config outputDirectory = get_climatology_op_directory(config) make_directories(outputDirectory) outputFileName = '{}/mocStreamfunction_years{:04d}-{:04d}.nc'.format( outputDirectory, self.startYear, self.endYear) if os.path.exists(outputFileName): return regionNames = config.getExpression(self.sectionName, 'regionNames') regionNames.append('Global') # Read in depth and bin latitudes try: restartFileName = self.runStreams.readpath('restart')[0] except ValueError: raise IOError('No MPAS-O restart file found: need at least ' 'one for MHT calcuation') with xr.open_dataset(restartFileName) as dsRestart: refBottomDepth = dsRestart.refBottomDepth.values nVertLevels = len(refBottomDepth) refLayerThickness = np.zeros(nVertLevels) refLayerThickness[0] = refBottomDepth[0] refLayerThickness[1:nVertLevels] = \ refBottomDepth[1:nVertLevels] - refBottomDepth[0:nVertLevels - 1] refZMid = refBottomDepth - 0.5 * refLayerThickness binBoundaryMocStreamfunction = None # first try timeSeriesStatsMonthly for bin boundaries, then try # mocStreamfunctionOutput stream as a backup option for streamName in ['timeSeriesStatsMonthlyOutput', 'mocStreamfunctionOutput']: try: inputFileName = self.historyStreams.readpath(streamName)[0] except ValueError: raise IOError('At least one file from stream {} is needed ' 'to compute MOC'.format(streamName)) with xr.open_dataset(inputFileName) as ds: if 'binBoundaryMocStreamfunction' in ds.data_vars: binBoundaryMocStreamfunction = \ ds.binBoundaryMocStreamfunction.values break if binBoundaryMocStreamfunction is None: raise ValueError('Could not find binBoundaryMocStreamfunction in ' 'either timeSeriesStatsMonthlyOutput or ' 'mocStreamfunctionOutput streams') binBoundaryMocStreamfunction = np.rad2deg(binBoundaryMocStreamfunction) # Compute and plot annual climatology of MOC streamfunction self.logger.info('\n Compute climatology of MOC streamfunction...') self.logger.info(' Load data...') climatologyFileName = self.mpasClimatologyTask.get_file_name( season='ANN') annualClimatology = xr.open_dataset(climatologyFileName) if 'Time' in annualClimatology.dims: annualClimatology = annualClimatology.isel(Time=0) # rename some variables for convenience annualClimatology = annualClimatology.rename( {'timeMonthly_avg_mocStreamvalLatAndDepth': 'avgMocStreamfunGlobal', 'timeMonthly_avg_mocStreamvalLatAndDepthRegion': 'avgMocStreamfunRegional'}) # Create dictionary for MOC climatology (NB: need this form # in order to convert it to xarray dataset later in the script) depth = refZMid lat = {} moc = {} for region in regionNames: self.logger.info(' Compute {} MOC...'.format(region)) if region == 'Global': mocTop = annualClimatology.avgMocStreamfunGlobal.values else: # hard-wire region=0 (Atlantic) for now indRegion = 0 mocVar = annualClimatology.avgMocStreamfunRegional mocTop = mocVar[indRegion, :, :].values # Store computed MOC to dictionary lat[region] = binBoundaryMocStreamfunction moc[region] = mocTop # Save to file self.logger.info(' Save global and regional MOC to file...') ncFile = netCDF4.Dataset(outputFileName, mode='w') # create dimensions ncFile.createDimension('nz', nVertLevels) for region in regionNames: latBins = lat[region] mocTop = moc[region] ncFile.createDimension('nx{}'.format(region), len(latBins)) # create variables x = ncFile.createVariable('lat{}'.format(region), 'f4', ('nx{}'.format(region),)) x.description = 'latitude bins for MOC {}'\ ' streamfunction'.format(region) x.units = 'degrees (-90 to 90)' y = ncFile.createVariable('moc{}'.format(region), 'f4', ('nz', 'nx{}'.format(region))) y.description = 'MOC {} streamfunction, annual'\ ' climatology'.format(region) y.units = 'Sv (10^6 m^3/s)' # save variables x[:] = latBins y[:, :] = mocTop depthVar = ncFile.createVariable('depth', 'f4', ('nz',)) depthVar.description = 'depth' depthVar.units = 'meters' depthVar[:] = depth ncFile.close() # }}} def _compute_moc_climo_postprocess(self): # {{{ """compute mean MOC streamfunction as a post-process""" config = self.config outputDirectory = get_climatology_op_directory(config) make_directories(outputDirectory) outputFileName = '{}/mocStreamfunction_years{:04d}-{:04d}.nc'.format( outputDirectory, self.startYear, self.endYear) if os.path.exists(outputFileName): return dvEdge, areaCell, refBottomDepth, latCell, nVertLevels, \ refTopDepth, refLayerThickness = _load_mesh(self.runStreams) regionNames = config.getExpression(self.sectionName, 'regionNames') # Load basin region related variables and save them to dictionary mpasMeshName = config.get('input', 'mpasMeshName') dictRegion = _build_region_mask_dict(config, regionNames, mpasMeshName, self.logger) # Add Global regionCellMask=1 everywhere to make the algorithm # for the global moc similar to that of the regional moc dictRegion['Global'] = { 'cellMask': np.ones(np.size(latCell))} regionNames.append('Global') # Compute and plot annual climatology of MOC streamfunction self.logger.info('\n Compute post-processed climatological of MOC ' 'streamfunction...') self.logger.info(' Load data...') climatologyFileName = self.mpasClimatologyTask.get_file_name( season='ANN') annualClimatology = xr.open_dataset(climatologyFileName) if 'Time' in annualClimatology.dims: annualClimatology = annualClimatology.isel(Time=0) if self.includeBolus: annualClimatology['avgNormalVelocity'] = \ annualClimatology['timeMonthly_avg_normalVelocity'] + \ annualClimatology['timeMonthly_avg_normalGMBolusVelocity'] annualClimatology['avgVertVelocityTop'] = \ annualClimatology['timeMonthly_avg_vertVelocityTop'] + \ annualClimatology['timeMonthly_avg_vertGMBolusVelocityTop'] else: # rename some variables for convenience annualClimatology = annualClimatology.rename( {'timeMonthly_avg_normalVelocity': 'avgNormalVelocity', 'timeMonthly_avg_vertVelocityTop': 'avgVertVelocityTop'}) # Convert to numpy arrays # (can result in a memory error for large array size) horizontalVel = annualClimatology.avgNormalVelocity.values verticalVel = annualClimatology.avgVertVelocityTop.values velArea = verticalVel * areaCell[:, np.newaxis] # Create dictionary for MOC climatology (NB: need this form # in order to convert it to xarray dataset later in the script) depth = refTopDepth lat = {} moc = {} for region in regionNames: self.logger.info(' Compute {} MOC...'.format(region)) self.logger.info(' Compute transport through region ' 'southern transect...') if region == 'Global': transportZ = np.zeros(nVertLevels) else: maxEdgesInTransect = \ dictRegion[region]['maxEdgesInTransect'] transectEdgeGlobalIDs = \ dictRegion[region]['transectEdgeGlobalIDs'] transectEdgeMaskSigns = \ dictRegion[region]['transectEdgeMaskSigns'] transportZ = _compute_transport(maxEdgesInTransect, transectEdgeGlobalIDs, transectEdgeMaskSigns, nVertLevels, dvEdge, refLayerThickness, horizontalVel) regionCellMask = dictRegion[region]['cellMask'] latBinSize = \ config.getfloat('streamfunctionMOC{}'.format(region), 'latBinSize') if region == 'Global': latBins = np.arange(-90.0, 90.1, latBinSize) else: indRegion = dictRegion[region]['indices'] latBins = latCell[indRegion] latBins = np.arange(np.amin(latBins), np.amax(latBins) + latBinSize, latBinSize) mocTop = _compute_moc(latBins, nVertLevels, latCell, regionCellMask, transportZ, velArea) # Store computed MOC to dictionary lat[region] = latBins moc[region] = mocTop # Save to file self.logger.info(' Save global and regional MOC to file...') ncFile = netCDF4.Dataset(outputFileName, mode='w') # create dimensions ncFile.createDimension('nz', len(refTopDepth)) for region in regionNames: latBins = lat[region] mocTop = moc[region] ncFile.createDimension('nx{}'.format(region), len(latBins)) # create variables x = ncFile.createVariable('lat{}'.format(region), 'f4', ('nx{}'.format(region),)) x.description = 'latitude bins for MOC {}'\ ' streamfunction'.format(region) x.units = 'degrees (-90 to 90)' y = ncFile.createVariable('moc{}'.format(region), 'f4', ('nz', 'nx{}'.format(region))) y.description = 'MOC {} streamfunction, annual'\ ' climatology'.format(region) y.units = 'Sv (10^6 m^3/s)' # save variables x[:] = latBins y[:, :] = mocTop depthVar = ncFile.createVariable('depth', 'f4', ('nz',)) depthVar.description = 'depth' depthVar.units = 'meters' depthVar[:] = depth ncFile.close() # }}} # }}} class PlotMOCClimatologySubtask(AnalysisTask): # {{{ """ Computation of a climatology of the model meridional overturning circulation. """ # Authors # ------- # Milena Veneziani, Mark Petersen, Phillip Wolfram, Xylar Asay-Davis def __init__(self, parentTask, controlConfig): # {{{ """ Construct the analysis task. Parameters ---------- parentTask : ``StreamfunctionMOC`` The main task of which this is a subtask controlConfig : ``MpasAnalysisConfigParser``, optional Configuration options for a control run (if any) """ # Authors # ------- # Xylar Asay-Davis # first, call the constructor from the base class (AnalysisTask) super(PlotMOCClimatologySubtask, self).__init__( config=parentTask.config, taskName=parentTask.taskName, componentName=parentTask.componentName, tags=parentTask.tags, subtaskName='plotMOCClimatology') parentTask.add_subtask(self) self.controlConfig = controlConfig # }}} def setup_and_check(self): # {{{ """ Perform steps to set up the analysis and check for errors in the setup. Raises ------ ValueError if timeSeriesStatsMonthly is not enabled in the MPAS run """ # Authors # ------- # Xylar Asay-Davis # first, call setup_and_check from the base class (AnalysisTask), # which will perform some common setup, including storing: # self.runDirectory , self.historyDirectory, self.plotsDirectory, # self.namelist, self.runStreams, self.historyStreams, # self.calendar super(PlotMOCClimatologySubtask, self).setup_and_check() config = self.config self.startYear = config.getint('climatology', 'startYear') self.endYear = config.getint('climatology', 'endYear') self.sectionName = 'streamfunctionMOC' self.xmlFileNames = [] self.filePrefixes = {} mainRunName = config.get('runs', 'mainRunName') self.regionNames = ['Global'] + config.getExpression(self.sectionName, 'regionNames') for region in self.regionNames: filePrefix = 'moc{}_{}_years{:04d}-{:04d}'.format( region, mainRunName, self.startYear, self.endYear) self.xmlFileNames.append('{}/{}.xml'.format(self.plotsDirectory, filePrefix)) self.filePrefixes[region] = filePrefix # }}} def run_task(self): # {{{ """ Plot the MOC climatology """ # Authors # ------- # Milena Veneziani, Mark Petersen, Phillip J. Wolfram, Xylar Asay-Davis self.logger.info("\nPlotting streamfunction of Meridional Overturning " "Circulation (MOC)...") config = self.config depth, lat, moc = self._load_moc(config) if self.controlConfig is None: refTitle = None diffTitle = None else: refDepth, refLat, refMOC = self._load_moc(self.controlConfig) refTitle = self.controlConfig.get('runs', 'mainRunName') diffTitle = 'Main - Control' # **** Plot MOC **** # Define plotting variables mainRunName = config.get('runs', 'mainRunName') movingAveragePointsClimatological = config.getint( self.sectionName, 'movingAveragePointsClimatological') colorbarLabel = '[Sv]' xLabel = 'latitude [deg]' yLabel = 'depth [m]' for region in self.regionNames: self.logger.info(' Plot climatological {} MOC...'.format(region)) title = '{} MOC (ANN, years {:04d}-{:04d})'.format( region, self.startYear, self.endYear) filePrefix = self.filePrefixes[region] outFileName = '{}/{}.png'.format(self.plotsDirectory, filePrefix) x = lat[region] z = depth regionMOC = moc[region] # Subset lat range minLat = config.getExpression('streamfunctionMOC{}'.format(region), 'latBinMin') maxLat = config.getExpression('streamfunctionMOC{}'.format(region), 'latBinMax') indLat = np.logical_and(x >= minLat, x <= maxLat) x = x[indLat] regionMOC = regionMOC[:, indLat] if self.controlConfig is None: refRegionMOC = None diff = None else: # the coords of the ref MOC won't necessarily match this MOC # so we need to interpolate nz, nx = regionMOC.shape refNz, refNx = refMOC[region].shape temp = np.zeros((refNz, nx)) for zIndex in range(refNz): temp[zIndex, :] = np.interp( x, refLat[region], refMOC[region][zIndex, :], left=np.nan, right=np.nan) refRegionMOC = np.zeros((nz, nx)) for xIndex in range(nx): refRegionMOC[:, xIndex] = np.interp( depth, refDepth, temp[:, xIndex], left=np.nan, right=np.nan) diff = regionMOC - refRegionMOC plot_vertical_section_comparison( config, x, z, regionMOC, refRegionMOC, diff, colorMapSectionName='streamfunctionMOC{}'.format(region), colorbarLabel=colorbarLabel, title=title, modelTitle=mainRunName, refTitle=refTitle, diffTitle=diffTitle, xlabel=xLabel, ylabel=yLabel, movingAveragePoints=movingAveragePointsClimatological) savefig(outFileName) caption = '{} Meridional Overturning Streamfunction'.format(region) write_image_xml( config=config, filePrefix=filePrefix, componentName='Ocean', componentSubdirectory='ocean', galleryGroup='Meridional Overturning Streamfunction', groupLink='moc', thumbnailDescription=region, imageDescription=caption, imageCaption=caption) # }}} def _load_moc(self, config): # {{{ """compute mean MOC streamfunction from analysis member""" startYear = config.getint('climatology', 'startYear') endYear = config.getint('climatology', 'endYear') outputDirectory = get_climatology_op_directory(config) make_directories(outputDirectory) inputFileName = '{}/mocStreamfunction_years{:04d}-{:04d}.nc'.format( outputDirectory, startYear, endYear) # Read from file ncFile = netCDF4.Dataset(inputFileName, mode='r') depth = ncFile.variables['depth'][:] lat = {} moc = {} for region in self.regionNames: lat[region] = ncFile.variables['lat{}'.format(region)][:] moc[region] = \ ncFile.variables['moc{}'.format(region)][:, :] ncFile.close() return depth, lat, moc # }}} # }}} class ComputeMOCTimeSeriesSubtask(AnalysisTask): # {{{ """ Computation of a time series of max Atlantic MOC at 26.5N. """ # Authors # ------- # Milena Veneziani, Mark Petersen, Phillip Wolfram, Xylar Asay-Davis def __init__(self, parentTask, startYear, endYear): # {{{ """ Construct the analysis task. Parameters ---------- parentTask : ``StreamfunctionMOC`` The main task of which this is a subtask """ # Authors # ------- # Xylar Asay-Davis # first, call the constructor from the base class (AnalysisTask) super(ComputeMOCTimeSeriesSubtask, self).__init__( config=parentTask.config, taskName=parentTask.taskName, componentName=parentTask.componentName, tags=parentTask.tags, subtaskName='computeMOCTimeSeries_{:04d}-{:04d}'.format( startYear, endYear)) parentTask.add_subtask(self) self.startYear = startYear self.endYear = endYear # }}} def setup_and_check(self): # {{{ """ Perform steps to set up the analysis and check for errors in the setup. Raises ------ ValueError if timeSeriesStatsMonthly is not enabled in the MPAS run """ # Authors # ------- # Xylar Asay-Davis # first, call setup_and_check from the base class (AnalysisTask), # which will perform some common setup, including storing: # self.runDirectory , self.historyDirectory, self.plotsDirectory, # self.namelist, self.runStreams, self.historyStreams, # self.calendar super(ComputeMOCTimeSeriesSubtask, self).setup_and_check() config = self.config self.mocAnalysisMemberEnabled = self.check_analysis_enabled( analysisOptionName='config_am_mocstreamfunction_enable', raiseException=False) self.sectionName = 'streamfunctionMOC' self.usePostprocessing = config.getExpression( self.sectionName, 'usePostprocessingScript') if not self.usePostprocessing and self.mocAnalysisMemberEnabled: self.variableList = \ ['timeMonthly_avg_mocStreamvalLatAndDepth', 'timeMonthly_avg_mocStreamvalLatAndDepthRegion'] else: self.variableList = ['timeMonthly_avg_normalVelocity', 'timeMonthly_avg_vertVelocityTop'] # Add the bolus velocity if GM is enabled try: # the new name self.includeBolus = self.namelist.getbool('config_use_gm') except KeyError: # the old name self.includeBolus = self.namelist.getbool( 'config_use_standardgm') if self.includeBolus: self.variableList.extend( ['timeMonthly_avg_normalGMBolusVelocity', 'timeMonthly_avg_vertGMBolusVelocityTop']) # }}} def run_task(self): # {{{ """ Process MOC analysis member data if available, or compute MOC at post-processing if not. """ # Authors # ------- # Milena Veneziani, Mark Petersen, Phillip J. Wolfram, Xylar Asay-Davis self.logger.info("\nCompute time series of Meridional Overturning " "Circulation (MOC)...") self.startDate = '{:04d}-01-01_00:00:00'.format(self.startYear) self.endDate = '{:04d}-12-31_23:59:59'.format(self.endYear) # **** Compute MOC **** if not self.usePostprocessing and self.mocAnalysisMemberEnabled: self._compute_moc_time_series_analysismember() else: self._compute_moc_time_series_postprocess() # }}} def _compute_moc_time_series_analysismember(self): # {{{ """compute MOC time series from analysis member""" # Compute and plot time series of Atlantic MOC at 26.5N (RAPID array) self.logger.info('\n Compute Atlantic MOC time series from analysis ' 'member...') self.logger.info(' Load data...') outputDirectory = '{}/moc/'.format( build_config_full_path(self.config, 'output', 'timeseriesSubdirectory')) try: os.makedirs(outputDirectory) except OSError: pass outputFileName = '{}/mocTimeSeries_{:04d}-{:04d}.nc'.format( outputDirectory, self.startYear, self.endYear) streamName = 'timeSeriesStatsMonthlyOutput' # Get bin latitudes and index of 26.5N binBoundaryMocStreamfunction = None # first try timeSeriesStatsMonthly for bin boundaries, then try # mocStreamfunctionOutput stream as a backup option for streamName in ['timeSeriesStatsMonthlyOutput', 'mocStreamfunctionOutput']: try: inputFileName = self.historyStreams.readpath(streamName)[0] except ValueError: raise IOError('At least one file from stream {} is needed ' 'to compute MOC'.format(streamName)) with xr.open_dataset(inputFileName) as ds: if 'binBoundaryMocStreamfunction' in ds.data_vars: binBoundaryMocStreamfunction = \ ds.binBoundaryMocStreamfunction.values break if binBoundaryMocStreamfunction is None: raise ValueError('Could not find binBoundaryMocStreamfunction in ' 'either timeSeriesStatsMonthlyOutput or ' 'mocStreamfunctionOutput streams') binBoundaryMocStreamfunction = np.rad2deg(binBoundaryMocStreamfunction) dLat = binBoundaryMocStreamfunction - 26.5 indlat26 = np.where(np.abs(dLat) == np.amin(np.abs(dLat))) inputFiles = sorted(self.historyStreams.readpath( streamName, startDate=self.startDate, endDate=self.endDate, calendar=self.calendar)) years, months = get_files_year_month(inputFiles, self.historyStreams, 'timeSeriesStatsMonthlyOutput') mocRegion = np.zeros(len(inputFiles)) times = np.zeros(len(inputFiles)) computed = np.zeros(len(inputFiles), bool) continueOutput = os.path.exists(outputFileName) if continueOutput: self.logger.info(' Read in previously computed MOC time series') with open_mpas_dataset(fileName=outputFileName, calendar=self.calendar, timeVariableNames=None, variableList=['mocAtlantic26'], startDate=self.startDate, endDate=self.endDate) as dsMOCIn: dsMOCIn.load() # first, copy all computed data for inIndex in range(dsMOCIn.dims['Time']): mask = np.logical_and( dsMOCIn.year[inIndex].values == years, dsMOCIn.month[inIndex].values == months) outIndex = np.where(mask)[0][0] mocRegion[outIndex] = dsMOCIn.mocAtlantic26[inIndex] times[outIndex] = dsMOCIn.Time[inIndex] computed[outIndex] = True if np.all(computed): # no need to waste time writing out the data set again return dsMOCIn for timeIndex, fileName in enumerate(inputFiles): if computed[timeIndex]: continue dsLocal = open_mpas_dataset( fileName=fileName, calendar=self.calendar, variableList=self.variableList, startDate=self.startDate, endDate=self.endDate) dsLocal = dsLocal.isel(Time=0) time = dsLocal.Time.values times[timeIndex] = time date = days_to_datetime(time, calendar=self.calendar) self.logger.info(' date: {:04d}-{:02d}'.format(date.year, date.month)) # hard-wire region=0 (Atlantic) for now indRegion = 0 mocVar = dsLocal.timeMonthly_avg_mocStreamvalLatAndDepthRegion mocTop = mocVar[indRegion, :, :].values mocRegion[timeIndex] = np.amax(mocTop[:, indlat26]) description = 'Max MOC Atlantic streamfunction nearest to RAPID ' \ 'Array latitude (26.5N)' dictonary = {'dims': ['Time'], 'coords': {'Time': {'dims': ('Time'), 'data': times, 'attrs': {'units': 'days since 0001-01-01'}}, 'year': {'dims': ('Time'), 'data': years, 'attrs': {'units': 'year'}}, 'month': {'dims': ('Time'), 'data': months, 'attrs': {'units': 'month'}}}, 'data_vars': {'mocAtlantic26': {'dims': ('Time'), 'data': mocRegion, 'attrs': {'units': 'Sv (10^6 m^3/s)', 'description': description}}}} dsMOCTimeSeries = xr.Dataset.from_dict(dictonary) write_netcdf(dsMOCTimeSeries, outputFileName) # }}} def _compute_moc_time_series_postprocess(self): # {{{ """compute MOC time series as a post-process""" config = self.config # Compute and plot time series of Atlantic MOC at 26.5N (RAPID array) self.logger.info('\n Compute and/or plot post-processed Atlantic MOC ' 'time series...') self.logger.info(' Load data...') outputDirectory = '{}/moc/'.format( build_config_full_path(self.config, 'output', 'timeseriesSubdirectory')) try: os.makedirs(outputDirectory) except OSError: pass outputFileName = '{}/mocTimeSeries_{:04d}-{:04d}.nc'.format( outputDirectory, self.startYear, self.endYear) dvEdge, areaCell, refBottomDepth, latCell, nVertLevels, \ refTopDepth, refLayerThickness = _load_mesh(self.runStreams) mpasMeshName = config.get('input', 'mpasMeshName') dictRegion = _build_region_mask_dict(config, ['Atlantic'], mpasMeshName, self.logger) dictRegion = dictRegion['Atlantic'] latBinSize = config.getfloat('streamfunctionMOCAtlantic', 'latBinSize') indRegion = dictRegion['indices'] latBins = latCell[indRegion] latBins = np.arange(np.amin(latBins), np.amax(latBins) + latBinSize, latBinSize) latAtlantic = latBins dLat = latAtlantic - 26.5 indlat26 = np.where(np.abs(dLat) == np.amin(np.abs(dLat))) maxEdgesInTransect = dictRegion['maxEdgesInTransect'] transectEdgeGlobalIDs = dictRegion['transectEdgeGlobalIDs'] transectEdgeMaskSigns = dictRegion['transectEdgeMaskSigns'] regionCellMask = dictRegion['cellMask'] streamName = 'timeSeriesStatsMonthlyOutput' inputFiles = sorted(self.historyStreams.readpath( streamName, startDate=self.startDate, endDate=self.endDate, calendar=self.calendar)) years, months = get_files_year_month(inputFiles, self.historyStreams, 'timeSeriesStatsMonthlyOutput') mocRegion = np.zeros(len(inputFiles)) times = np.zeros(len(inputFiles)) computed = np.zeros(len(inputFiles), bool) continueOutput = os.path.exists(outputFileName) if continueOutput: self.logger.info(' Read in previously computed MOC time series') with open_mpas_dataset(fileName=outputFileName, calendar=self.calendar, timeVariableNames=None, variableList=['mocAtlantic26'], startDate=self.startDate, endDate=self.endDate) as dsMOCIn: dsMOCIn.load() # first, copy all computed data for inIndex in range(dsMOCIn.dims['Time']): mask = np.logical_and( dsMOCIn.year[inIndex].values == years, dsMOCIn.month[inIndex].values == months) outIndex = np.where(mask)[0][0] mocRegion[outIndex] = dsMOCIn.mocAtlantic26[inIndex] times[outIndex] = dsMOCIn.Time[inIndex] computed[outIndex] = True if np.all(computed): # no need to waste time writing out the data set again return dsMOCIn for timeIndex, fileName in enumerate(inputFiles): if computed[timeIndex]: continue dsLocal = open_mpas_dataset( fileName=fileName, calendar=self.calendar, variableList=self.variableList, startDate=self.startDate, endDate=self.endDate) dsLocal = dsLocal.isel(Time=0) time = dsLocal.Time.values times[timeIndex] = time date = days_to_datetime(time, calendar=self.calendar) self.logger.info(' date: {:04d}-{:02d}'.format(date.year, date.month)) if self.includeBolus: dsLocal['avgNormalVelocity'] = \ dsLocal['timeMonthly_avg_normalVelocity'] + \ dsLocal['timeMonthly_avg_normalGMBolusVelocity'] dsLocal['avgVertVelocityTop'] = \ dsLocal['timeMonthly_avg_vertVelocityTop'] + \ dsLocal['timeMonthly_avg_vertGMBolusVelocityTop'] else: # rename some variables for convenience dsLocal = dsLocal.rename( {'timeMonthly_avg_normalVelocity': 'avgNormalVelocity', 'timeMonthly_avg_vertVelocityTop': 'avgVertVelocityTop'}) horizontalVel = dsLocal.avgNormalVelocity.values verticalVel = dsLocal.avgVertVelocityTop.values velArea = verticalVel * areaCell[:, np.newaxis] transportZ = _compute_transport(maxEdgesInTransect, transectEdgeGlobalIDs, transectEdgeMaskSigns, nVertLevels, dvEdge, refLayerThickness, horizontalVel) mocTop = _compute_moc(latAtlantic, nVertLevels, latCell, regionCellMask, transportZ, velArea) mocRegion[timeIndex] = np.amax(mocTop[:, indlat26]) description = 'Max MOC Atlantic streamfunction nearest to RAPID ' \ 'Array latitude (26.5N)' dictonary = {'dims': ['Time'], 'coords': {'Time': {'dims': ('Time'), 'data': times, 'attrs': {'units': 'days since 0001-01-01'}}, 'year': {'dims': ('Time'), 'data': years, 'attrs': {'units': 'year'}}, 'month': {'dims': ('Time'), 'data': months, 'attrs': {'units': 'month'}}}, 'data_vars': {'mocAtlantic26': {'dims': ('Time'), 'data': mocRegion, 'attrs': {'units': 'Sv (10^6 m^3/s)', 'description': description}}}} dsMOCTimeSeries = xr.Dataset.from_dict(dictonary) write_netcdf(dsMOCTimeSeries, outputFileName) # }}} # }}} class CombineMOCTimeSeriesSubtask(AnalysisTask): # {{{ """ Combine individual time series of max Atlantic MOC at 26.5N into a single data set """ # Authors # ------- # Xylar Asay-Davis def __init__(self, parentTask, startYears, endYears): # {{{ """ Construct the analysis task. Parameters ---------- parentTask : ``StreamfunctionMOC`` The main task of which this is a subtask controlConfig : ``MpasAnalysisConfigParser``, optional Configuration options for a control run (if any) """ # Authors # ------- # Xylar Asay-Davis # first, call the constructor from the base class (AnalysisTask) super(CombineMOCTimeSeriesSubtask, self).__init__( config=parentTask.config, taskName=parentTask.taskName, componentName=parentTask.componentName, tags=parentTask.tags, subtaskName='combineMOCTimeSeries') parentTask.add_subtask(self) self.startYears = startYears self.endYears = endYears # }}} def run_task(self): # {{{ """ Plot the MOC time series """ # Authors # ------- # Xylar Asay-Davis outputDirectory = '{}/moc/'.format( build_config_full_path(self.config, 'output', 'timeseriesSubdirectory')) try: os.makedirs(outputDirectory) except OSError: pass outputFileNames = [] for startYear, endYear in zip(self.startYears, self.endYears): outputFileName = '{}/mocTimeSeries_{:04d}-{:04d}.nc'.format( outputDirectory, startYear, endYear) outputFileNames.append(outputFileName) outputFileName = '{}/mocTimeSeries_{:04d}-{:04d}.nc'.format( outputDirectory, self.startYears[0], self.endYears[-1]) if outputFileName in outputFileNames: # don't try to write to read from and write to the same file return ds = xr.open_mfdataset(outputFileNames, concat_dim='Time', combine='nested', decode_times=False) ds.load() write_netcdf(ds, outputFileName) # }}} # }}} class PlotMOCTimeSeriesSubtask(AnalysisTask): # {{{ """ Plots a time series of max Atlantic MOC at 26.5N. """ # Authors # ------- # Milena Veneziani, Mark Petersen, Phillip Wolfram, Xylar Asay-Davis def __init__(self, parentTask, controlConfig): # {{{ """ Construct the analysis task. Parameters ---------- parentTask : ``StreamfunctionMOC`` The main task of which this is a subtask controlConfig : ``MpasAnalysisConfigParser``, optional Configuration options for a control run (if any) """ # Authors # ------- # Xylar Asay-Davis # first, call the constructor from the base class (AnalysisTask) super(PlotMOCTimeSeriesSubtask, self).__init__( config=parentTask.config, taskName=parentTask.taskName, componentName=parentTask.componentName, tags=parentTask.tags, subtaskName='plotMOCTimeSeries') parentTask.add_subtask(self) self.controlConfig = controlConfig # }}} def setup_and_check(self): # {{{ """ Perform steps to set up the analysis and check for errors in the setup. Raises ------ ValueError if timeSeriesStatsMonthly is not enabled in the MPAS run """ # Authors # ------- # Xylar Asay-Davis # first, call setup_and_check from the base class (AnalysisTask), # which will perform some common setup, including storing: # self.runDirectory , self.historyDirectory, self.plotsDirectory, # self.namelist, self.runStreams, self.historyStreams, # self.calendar super(PlotMOCTimeSeriesSubtask, self).setup_and_check() config = self.config self.sectionName = 'streamfunctionMOC' mainRunName = config.get('runs', 'mainRunName') filePrefix = 'mocTimeseries_{}'.format(mainRunName) self.xmlFileNames = ['{}/{}.xml'.format(self.plotsDirectory, filePrefix)] self.filePrefix = filePrefix # }}} def run_task(self): # {{{ """ Plot the MOC time series """ # Authors # ------- # Milena Veneziani, Mark Petersen, Phillip J. Wolfram, Xylar Asay-Davis self.logger.info("\nPlotting time series of Meridional Overturning " "Circulation (MOC)...") config = self.config dsMOCTimeSeries = self._load_moc(config) # **** Plot MOC **** # Define plotting variables mainRunName = config.get('runs', 'mainRunName') movingAveragePoints = config.getint(self.sectionName, 'movingAveragePoints') # Plot time series self.logger.info(' Plot time series of max Atlantic MOC at 26.5N...') xLabel = 'Time [years]' yLabel = '[Sv]' title = r'Max Atlantic MOC at $26.5\degree$N\n {}'.format(mainRunName) filePrefix = self.filePrefix outFileName = '{}/{}.png'.format(self.plotsDirectory, filePrefix) if config.has_option(self.taskName, 'firstYearXTicks'): firstYearXTicks = config.getint(self.taskName, 'firstYearXTicks') else: firstYearXTicks = None if config.has_option(self.taskName, 'yearStrideXTicks'): yearStrideXTicks = config.getint(self.taskName, 'yearStrideXTicks') else: yearStrideXTicks = None fields = [dsMOCTimeSeries.mocAtlantic26] lineColors = ['k'] lineWidths = [2] legendText = [mainRunName] if self.controlConfig is not None: dsRefMOC = self._load_moc(self.controlConfig) fields.append(dsRefMOC.mocAtlantic26) lineColors.append('r') lineWidths.append(2) controlRunName = self.controlConfig.get('runs', 'mainRunName') legendText.append(controlRunName) timeseries_analysis_plot(config, fields, calendar=self.calendar, title=title, xlabel=xLabel, ylabel=yLabel, movingAveragePoints=movingAveragePoints, lineColors=lineColors, lineWidths=lineWidths, legendText=legendText, firstYearXTicks=firstYearXTicks, yearStrideXTicks=yearStrideXTicks) savefig(outFileName) caption = u'Time Series of maximum Meridional Overturning ' \ u'Circulation at 26.5°N' write_image_xml( config=config, filePrefix=filePrefix, componentName='Ocean', componentSubdirectory='ocean', galleryGroup='Meridional Overturning Streamfunction', groupLink='moc', thumbnailDescription='Time Series', imageDescription=caption, imageCaption=caption) # }}} def _load_moc(self, config): # {{{ """compute mean MOC streamfunction from analysis member""" outputDirectory = build_config_full_path(config, 'output', 'timeseriesSubdirectory') startYear = config.getint('timeSeries', 'startYear') endYear = config.getint('timeSeries', 'endYear') inputFileName = '{}/moc/mocTimeSeries_{:04d}-{:04d}.nc'.format( outputDirectory, startYear, endYear) dsMOCTimeSeries = xr.open_dataset(inputFileName, decode_times=False) return dsMOCTimeSeries # }}} # }}} def _load_mesh(runStreams): # {{{ # Load mesh related variables try: restartFile = runStreams.readpath('restart')[0] except ValueError: raise IOError('No MPAS-O restart file found: need at least one ' 'restart file for MOC calculation') ncFile = netCDF4.Dataset(restartFile, mode='r') dvEdge = ncFile.variables['dvEdge'][:] areaCell = ncFile.variables['areaCell'][:] refBottomDepth = ncFile.variables['refBottomDepth'][:] latCell = np.rad2deg(ncFile.variables['latCell'][:]) ncFile.close() nVertLevels = len(refBottomDepth) refTopDepth = np.zeros(nVertLevels + 1) refTopDepth[1:nVertLevels + 1] = refBottomDepth[0:nVertLevels] refLayerThickness = np.zeros(nVertLevels) refLayerThickness[0] = refBottomDepth[0] refLayerThickness[1:nVertLevels] = \ (refBottomDepth[1:nVertLevels] - refBottomDepth[0:nVertLevels - 1]) return dvEdge, areaCell, refBottomDepth, latCell, nVertLevels, \ refTopDepth, refLayerThickness # }}} def _build_region_mask_dict(config, regionNames, mpasMeshName, logger): # {{{ regionMaskFile = get_region_mask( config, '{}_moc_masks_and_transects.nc'.format(mpasMeshName)) if not os.path.exists(regionMaskFile): raise IOError('Regional masking file {} for MOC calculation ' 'does not exist'.format(regionMaskFile)) iRegion = 0 dictRegion = {} for region in regionNames: logger.info('\n Reading region and transect mask for ' '{}...'.format(region)) ncFileRegional = netCDF4.Dataset(regionMaskFile, mode='r') maxEdgesInTransect = \ ncFileRegional.dimensions['maxEdgesInTransect'].size transectEdgeMaskSigns = \ ncFileRegional.variables['transectEdgeMaskSigns'][:, iRegion] transectEdgeGlobalIDs = \ ncFileRegional.variables['transectEdgeGlobalIDs'][iRegion, :] regionCellMask = \ ncFileRegional.variables['regionCellMasks'][:, iRegion] ncFileRegional.close() iRegion += 1 indRegion = np.where(regionCellMask == 1) dictRegion[region] = { 'indices': indRegion, 'cellMask': regionCellMask, 'maxEdgesInTransect': maxEdgesInTransect, 'transectEdgeMaskSigns': transectEdgeMaskSigns, 'transectEdgeGlobalIDs': transectEdgeGlobalIDs} return dictRegion # }}} def _compute_transport(maxEdgesInTransect, transectEdgeGlobalIDs, transectEdgeMaskSigns, nz, dvEdge, refLayerThickness, horizontalVel): # {{{ """compute mass transport across southern transect of ocean basin""" transportZEdge = np.zeros([nz, maxEdgesInTransect]) for i in range(maxEdgesInTransect): if transectEdgeGlobalIDs[i] == 0: break # subtract 1 because of python 0-indexing iEdge = transectEdgeGlobalIDs[i] - 1 transportZEdge[:, i] = horizontalVel[iEdge, :] * \ transectEdgeMaskSigns[iEdge, np.newaxis] * \ dvEdge[iEdge, np.newaxis] * \ refLayerThickness[np.newaxis, :] transportZ = transportZEdge.sum(axis=1) return transportZ # }}} def _compute_moc(latBins, nz, latCell, regionCellMask, transportZ, velArea): # {{{ """compute meridionally integrated MOC streamfunction""" mocTop = np.zeros([np.size(latBins), nz + 1]) mocTop[0, range(1, nz + 1)] = transportZ.cumsum() for iLat in range(1, np.size(latBins)): indlat = np.logical_and(np.logical_and( regionCellMask == 1, latCell >= latBins[iLat - 1]), latCell < latBins[iLat]) mocTop[iLat, :] = mocTop[iLat - 1, :] + \ velArea[indlat, :].sum(axis=0) # convert m^3/s to Sverdrup mocTop = mocTop * m3ps_to_Sv mocTop = mocTop.T return mocTop # }}} # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python