diff --git a/mpas_analysis/ocean/climatology_map_antarctic_melt.py b/mpas_analysis/ocean/climatology_map_antarctic_melt.py index 4c02cfad3..f01ce8be8 100644 --- a/mpas_analysis/ocean/climatology_map_antarctic_melt.py +++ b/mpas_analysis/ocean/climatology_map_antarctic_melt.py @@ -76,7 +76,6 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, sectionName = self.taskName - mpasFieldName = 'timeMonthly_avg_landIceFreshwaterFlux' iselValues = None # read in what seasons we want to plot @@ -109,7 +108,7 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, mpasClimatologyTask=mpasClimatologyTask, parentTask=self, climatologyName=fieldName, - variableList=[mpasFieldName], + variableList=None, comparisonGridNames=comparisonGridNames, seasons=seasons, iselValues=iselValues) @@ -117,7 +116,7 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, if controlConfig is None: refTitleLabel = \ - 'Observations (Paolo et al, 2023)' + 'Observations (Paolo et al. 2023)' observationsDirectory = build_obs_path( config, 'ocean', 'meltSubdirectory') @@ -138,13 +137,10 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, f'{observationsDirectory}/Paolo/Paolo_2023_' \ f'iceshelf_melt_rates_1992-2017_v1.0_6000x6000km_{res:g}km_' \ f'Antarctic_stereo.20240220.nc' - refFieldName = 'meltRate' - outFileLabel = 'meltPaolo' - galleryName = 'Observations: Paolo et al. (2023)' remapObservationsSubtask = RemapObservedAntarcticMeltClimatology( parentTask=self, seasons=seasons, fileName=obsFileName, - outFilePrefix=refFieldName, + outFilePrefix='meltRate', comparisonGridNames=comparisonGridNames) self.add_subtask(remapObservationsSubtask) diffTitleLabel = 'Model - Observations' @@ -152,33 +148,117 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, else: remapObservationsSubtask = None controlRunName = controlConfig.get('runs', 'mainRunName') - galleryName = None refTitleLabel = f'Control: {controlRunName}' + diffTitleLabel = 'Main - Control' + + totalFluxVar = 'timeMonthly_avg_landIceFreshwaterFluxTotal' + landIceFluxVar = 'timeMonthly_avg_landIceFreshwaterFlux' + frazilFluxVar = 'timeMonthly_avg_frazilIceFreshwaterFlux' + mpasFieldName = totalFluxVar + + if controlConfig is None: + refFieldName = 'meltRate' + else: refFieldName = mpasFieldName - outFileLabel = 'melt' - diffTitleLabel = 'Main - Control' for comparisonGridName in comparisonGridNames: for season in seasons: # make a new subtask for this season and comparison grid - subtask = PlotClimatologyMapSubtask( + subtaskName = f'plot_total_melt_{season}_{comparisonGridName}' + subtask = PlotAntarcticMeltSubtask( self, season, comparisonGridName, remapClimatologySubtask, - remapObservationsSubtask, controlConfig=controlConfig) + remapObservationsSubtask, controlConfig=controlConfig, + subtaskName=subtaskName) subtask.set_plot_info( - outFileLabel=outFileLabel, - fieldNameInTitle='Melt Rate', + outFileLabel='antMeltTotal', + fieldNameInTitle='Total Melt Flux', mpasFieldName=mpasFieldName, refFieldName=refFieldName, refTitleLabel=refTitleLabel, diffTitleLabel=diffTitleLabel, - unitsLabel=r'm a$^{-1}$', - imageCaption='Antarctic Melt Rate', + unitsLabel=r'm a$^{-1}$ freshwater equiv.', + imageCaption='Antarctic Total Melt Flux', galleryGroup='Melt Rate', groupSubtitle=None, groupLink='antarctic_melt', - galleryName=galleryName) + galleryName='Total Melt Flux') + + self.add_subtask(subtask) + + mpasFieldName = landIceFluxVar + + if controlConfig is None: + refFieldName = 'meltRate' + else: + refFieldName = mpasFieldName + + for comparisonGridName in comparisonGridNames: + for season in seasons: + # make a new subtask for this season and comparison grid + subtaskName = \ + f'plot_interface_melt_{season}_{comparisonGridName}' + subtask = PlotAntarcticMeltSubtask( + self, season, comparisonGridName, remapClimatologySubtask, + remapObservationsSubtask, controlConfig=controlConfig, + subtaskName=subtaskName) + + # In PlotAntarcticMeltSubtask, we will remove the obs from + # these plots if totalFluxVar is present so we only compare one + # field with obs + + subtask.set_plot_info( + outFileLabel='antMeltInterface', + fieldNameInTitle='Melt Rate at Interface', + mpasFieldName=mpasFieldName, + refFieldName=refFieldName, + refTitleLabel=refTitleLabel, + diffTitleLabel=diffTitleLabel, + unitsLabel=r'm a$^{-1}$ freshwater equiv.', + imageCaption='Antarctic Melt Rate at Interface', + galleryGroup='Melt Rate', + groupSubtitle=None, + groupLink='antarctic_melt_int', + galleryName='Melt Rate at the Ice-ocean Interface') + + self.add_subtask(subtask) + + mpasFieldName = frazilFluxVar + + if controlConfig is None: + refTitleLabel = None + refFieldName = None + diffTitleLabel = None + + else: + controlRunName = controlConfig.get('runs', 'mainRunName') + refTitleLabel = f'Control: {controlRunName}' + refFieldName = mpasFieldName + diffTitleLabel = 'Main - Control' + + for comparisonGridName in comparisonGridNames: + for season in seasons: + # make a new subtask for this season and comparison grid + subtaskName = \ + f'plot_interface_frazil_{season}_{comparisonGridName}' + subtask = PlotAntarcticMeltSubtask( + self, season, comparisonGridName, remapClimatologySubtask, + controlConfig=controlConfig, subtaskName=subtaskName) + + subtask.set_plot_info( + outFileLabel='antFrazil', + fieldNameInTitle='Frazil Accretion Rate, negative upward', + mpasFieldName=mpasFieldName, + refFieldName=refFieldName, + refTitleLabel=refTitleLabel, + diffTitleLabel=diffTitleLabel, + unitsLabel=r'm a$^{-1}$ freshwater equiv.', + imageCaption='Antarctic Accretion Rate', + galleryGroup='Melt Rate', + groupSubtitle=None, + groupLink='antarctic_frazil_flux', + galleryName='Frazil Accretion Rate') self.add_subtask(subtask) @@ -212,11 +292,35 @@ class RemapMpasAntarcticMeltClimatology(RemapMpasClimatologySubtask): landIceMask : xarray.DataArray A mask indicating where there is land ice on the ocean grid (thus, where melt rates are valid) + + renameDict : dict + A dictionary use to rename variables in the climatology """ # Authors # ------- # Xylar Asay-Davis + def setup_and_check(self): + """ + Figure out which variable(s) to remap + """ + # Authors + # ------- + # Xylar Asay-Davis + + totalFluxVar = 'timeMonthly_avg_landIceFreshwaterFluxTotal' + landIceFluxVar = 'timeMonthly_avg_landIceFreshwaterFlux' + frazilFluxVar = 'timeMonthly_avg_frazilIceFreshwaterFlux' + + if totalFluxVar in self.mpasClimatologyTask.allVariables: + # include the total and constituent fluxes + self.variableList = [totalFluxVar, landIceFluxVar, frazilFluxVar] + else: + # we only have the old name without the frazil accretion rate + self.variableList = [landIceFluxVar] + + super().setup_and_check() + def run_task(self): """ Compute climatologies of melt rates from E3SM/MPAS output @@ -260,12 +364,14 @@ def customize_masked_climatology(self, climatology, season): # ------- # Xylar Asay-Davis - fieldName = self.variableList[0] + for fieldName in self.variableList: - # scale the field to m/yr from kg/m^2/s and mask out non-land-ice areas - climatology[fieldName] = \ - constants.sec_per_year / constants.rho_fw * \ - climatology[fieldName].where(self.landIceMask) + # scale the field to m/yr from kg/m^2/s and mask out non-land-ice + # areas + climatology[fieldName] = \ + constants.sec_per_year / constants.rho_fw * \ + climatology[fieldName].where(self.landIceMask) + climatology[fieldName].attrs['units'] = 'm yr^-1' return climatology @@ -422,7 +528,8 @@ def run_task(self): pool=ThreadPool(1)): # Load data: - inFileName = self.mpasClimatologyTask.get_file_name(self.season) + inFileName = \ + self.mpasClimatologyTask.get_file_name(self.season) mpasFieldName = 'timeMonthly_avg_landIceFreshwaterFlux' dsIn = xr.open_dataset(inFileName) freshwaterFlux = dsIn[mpasFieldName] @@ -445,7 +552,8 @@ def run_task(self): # select only those regions we want to plot dsRegionMask = dsRegionMask.isel(nRegions=regionIndices) - cellMasks = dsRegionMask.regionCellMasks.chunk({'nRegions': 10}) + cellMasks = \ + dsRegionMask.regionCellMasks.chunk({'nRegions': 10}) restartFileName = \ self.runStreams.readpath('restart')[0] @@ -548,3 +656,53 @@ def run_task(self): row[controlRunName] = \ f'{dsControl.totalMeltFlux[index].values}' writer.writerow(row) + + +class PlotAntarcticMeltSubtask(PlotClimatologyMapSubtask): + """ + A subtask for plotting antarctic melt fields if available + + Attributes + ---------- + doPlot : bool + Whether the required variable from the climatology is available so that + a plot should be generated + """ + # Authors + # ------- + # Xylar Asay-Davis + + def setup_and_check(self): + """ + Perform steps to set up the analysis and check for errors in the setup. + """ + allVariables = \ + self.remapMpasClimatologySubtask.mpasClimatologyTask.allVariables + + totalFluxVar = 'timeMonthly_avg_landIceFreshwaterFluxTotal' + landIceFluxVar = 'timeMonthly_avg_landIceFreshwaterFlux' + plotAll = (totalFluxVar in allVariables) + + if self.mpasFieldName == landIceFluxVar and plotAll and \ + self.controlConfig is None: + # need to remove obs because we only wnat to plot them vs the + # total flux + self.remapObsClimatologySubtask = None + self.refTitleLabel = None + self.refFieldName = None + self.diffTitleLabel = None + + self.doPlot = (self.mpasFieldName == landIceFluxVar or plotAll) + + if self.doPlot: + super().setup_and_check() + else: + # still need to call the base class's method + AnalysisTask.setup_and_check(self=self) + + def run_task(self): + """ + Plot the variable if available + """ + if self.doPlot: + super().run_task() diff --git a/mpas_analysis/ocean/time_series_antarctic_melt.py b/mpas_analysis/ocean/time_series_antarctic_melt.py index 4bffa8793..d509b5685 100644 --- a/mpas_analysis/ocean/time_series_antarctic_melt.py +++ b/mpas_analysis/ocean/time_series_antarctic_melt.py @@ -80,7 +80,8 @@ def __init__(self, config, mpasTimeSeriesTask, regionMasksTask, # nothing else to do return - masksSubtask = regionMasksTask.add_mask_subtask(regionGroup=regionGroup) + masksSubtask = \ + regionMasksTask.add_mask_subtask(regionGroup=regionGroup) self.iceShelfMasksFile = masksSubtask.geojsonFileName iceShelvesToPlot = masksSubtask.expand_region_names(iceShelvesToPlot) @@ -185,8 +186,7 @@ def __init__(self, parentTask, startYear, endYear, mpasTimeSeriesTask, self.endYear = endYear self.startDate = f'{self.startYear:04d}-01-01_00:00:00' self.endDate = f'{self.endYear:04d}-12-31_23:59:59' - self.variableList = \ - ['timeMonthly_avg_landIceFreshwaterFlux'] + self.variableList = None def setup_and_check(self): """ @@ -230,6 +230,13 @@ def setup_and_check(self): raise IOError('No MPAS-O restart file found: need at least one ' 'restart file for Antarctic melt calculations') + totalFluxVar = 'timeMonthly_avg_landIceFreshwaterFluxTotal' + landIceFluxVar = 'timeMonthly_avg_landIceFreshwaterFlux' + if totalFluxVar in self.mpasTimeSeriesTask.allVariables: + self.variableList = [totalFluxVar] + else: + self.variableList = [landIceFluxVar] + self.mpasTimeSeriesTask.add_variables(variableList=self.variableList) return @@ -309,17 +316,18 @@ def run_task(self): regionNames = decode_strings(dsRegionMask.regionNames) + fluxVar = self.variableList[0] + datasets = [] nTime = dsIn.sizes['Time'] for tIndex in range(nTime): self.logger.info(f' {tIndex + 1}/{nTime}') - freshwaterFlux = \ - dsIn.timeMonthly_avg_landIceFreshwaterFlux.isel(Time=tIndex) + freshwaterFlux = dsIn[fluxVar].isel(Time=tIndex) nRegions = dsRegionMask.sizes['nRegions'] meltRates = numpy.zeros((nRegions,)) - totalMeltFluxes = numpy.zeros((nRegions,)) + integratedMeltFluxes = numpy.zeros((nRegions,)) for regionIndex in range(nRegions): self.logger.info(f' {regionNames[regionIndex]}') @@ -327,7 +335,7 @@ def run_task(self): dsRegionMask.regionCellMasks.isel(nRegions=regionIndex) # convert from kg/s to kg/yr - totalMeltFlux = constants.sec_per_year * \ + integratedMeltFlux = constants.sec_per_year * \ (cellMask * areaCell * freshwaterFlux).sum(dim='nCells') totalArea = \ @@ -335,23 +343,23 @@ def run_task(self): # from kg/m^2/yr to m/yr meltRates[regionIndex] = ((1. / constants.rho_fw) * - (totalMeltFlux / totalArea)) + (integratedMeltFlux / totalArea)) # convert from kg/yr to GT/yr - totalMeltFlux /= constants.kg_per_GT - totalMeltFluxes[regionIndex] = totalMeltFlux + integratedMeltFlux /= constants.kg_per_GT + integratedMeltFluxes[regionIndex] = integratedMeltFlux dsOut = xarray.Dataset() dsOut.coords['Time'] = dsIn.Time.isel(Time=tIndex) - dsOut['totalMeltFlux'] = (('nRegions',), totalMeltFluxes) + dsOut['integratedMeltFlux'] = (('nRegions',), integratedMeltFluxes) dsOut['meltRates'] = (('nRegions',), meltRates) datasets.append(dsOut) dsOut = xarray.concat(objs=datasets, dim='Time') dsOut['regionNames'] = dsRegionMask.regionNames - dsOut.totalMeltFlux.attrs['units'] = 'GT a$^{-1}$' - dsOut.totalMeltFlux.attrs['description'] = \ - 'Total melt flux summed over each ice shelf or region' + dsOut.integratedMeltFlux.attrs['units'] = 'GT a$^{-1}$' + dsOut.integratedMeltFlux.attrs['description'] = \ + 'Integrated melt flux summed over each ice shelf or region' dsOut.meltRates.attrs['units'] = 'm a$^{-1}$' dsOut.meltRates.attrs['description'] = \ 'Melt rate averaged over each ice shelf or region' @@ -536,17 +544,17 @@ def run_task(self): fc.add_feature(feature) break - totalMeltFlux, meltRates = self._load_ice_shelf_fluxes(config) + integratedMeltFlux, meltRates = self._load_ice_shelf_fluxes(config) plotControl = self.controlConfig is not None if plotControl: controlRunName = self.controlConfig.get('runs', 'mainRunName') - refTotalMeltFlux, refMeltRates = \ + refintegratedMeltFlux, refMeltRates = \ self._load_ice_shelf_fluxes(self.controlConfig) else: controlRunName = None - refTotalMeltFlux = None + refintegratedMeltFlux = None refMeltRates = None # Load observations from multiple files and put in dictionary based @@ -595,7 +603,7 @@ def run_task(self): index = region_names.index(self.iceShelf) ds_shelf = ds_adusumilli.isel(nRegions=index) obsDict['Adusumilli et al. (2020)'] = { - 'meltFlux': ds_shelf.totalMeltFlux.values, + 'meltFlux': ds_shelf.integratedMeltFlux.values, 'meltFluxUncertainty': ds_shelf.meltFluxUncertainty.values, 'meltRate': ds_shelf.meanMeltRate.values, 'meltRateUncertainty': ds_shelf.meltRateUncertainty.values} @@ -671,7 +679,7 @@ def run_task(self): xLabel = 'Time (yr)' yLabel = 'Melt Flux (GT/yr)' - timeSeries = totalMeltFlux.isel(nRegions=self.regionIndex) + timeSeries = integratedMeltFlux.isel(nRegions=self.regionIndex) filePrefix = f'melt_flux_{suffix}' outFileName = f'{self.plotsDirectory}/{filePrefix}.png' @@ -681,7 +689,7 @@ def run_task(self): lineWidths = [2.5] legendText = [mainRunName] if plotControl: - fields.append(refTotalMeltFlux.isel(nRegions=self.regionIndex)) + fields.append(refintegratedMeltFlux.isel(nRegions=self.regionIndex)) lineColors.append(config.get('timeSeries', 'controlColor')) lineWidths.append(1.2) legendText.append(controlRunName) @@ -734,7 +742,7 @@ def run_task(self): savefig(outFileName, config) - caption = f'Running Mean of Total Melt Flux under Ice Shelves in ' \ + caption = f'Running Mean of Integrated Melt Flux under Ice Shelves in ' \ f'the {title} Region' write_image_xml( config=config, @@ -743,13 +751,13 @@ def run_task(self): componentSubdirectory='ocean', galleryGroup='Antarctic Melt Time Series', groupLink='antmelttime', - gallery='Total Melt Flux', + gallery='Integrated Melt Flux', thumbnailDescription=title, imageDescription=caption, imageCaption=caption) xLabel = 'Time (yr)' - yLabel = 'Melt Rate (m/yr)' + yLabel = 'Melt Rate (m/yr) freshwater equiv.' timeSeries = meltRates.isel(nRegions=self.regionIndex) @@ -816,8 +824,8 @@ def run_task(self): @staticmethod def _load_ice_shelf_fluxes(config): """ - Reads melt flux time series and computes regional total melt flux and - mean melt rate. + Reads melt flux time series and computes regional integrated melt flux + and mean melt rate. """ # Authors # ------- @@ -835,4 +843,4 @@ def _load_ice_shelf_fluxes(config): f'{startYear:04d}-{endYear:04d}.nc' dsOut = xarray.open_dataset(outFileName) - return dsOut.totalMeltFlux, dsOut.meltRates + return dsOut.integratedMeltFlux, dsOut.meltRates