# coding=utf-8
"""Compute the total ocean heat content"""
import shutil
import numpy as np
from earthdiagnostics import cdftools
from earthdiagnostics.box import Box
from earthdiagnostics.constants import Basins
from earthdiagnostics.diagnostic import Diagnostic, DiagnosticBasinOption, DiagnosticIntOption
from earthdiagnostics.modelingrealm import ModelingRealms
from earthdiagnostics.utils import Utils, TempFile
[docs]class HeatContent(Diagnostic):
"""
Compute the total ocean heat content
:original author: Virginie Guemas <virginie.guemas@bsc.es>
:contributor: Javier Vegas-Regidor<javier.vegas@bsc.es>
:created: May 2012
:last modified: June 2016
:param data_manager: data management object
:type data_manager: DataManager
:param startdate: startdate
:type startdate: str
:param member: member number
:type member: int
:param chunk: chunk's number
:type chunk: int
:param mixed_layer: If 1, restricts calculation to the mixed layer, if -1 exclude it. If 0, no effect
:type mixed_layer: int
:param box: box to use for the average
:type box: Box
"""
alias = 'ohc'
"Diagnostic alias for the configuration file"
def __init__(self, data_manager, startdate, member, chunk, basin, mixed_layer, box, min_level, max_level):
Diagnostic.__init__(self, data_manager)
self.startdate = startdate
self.member = member
self.chunk = chunk
self.basin = basin
self.mxloption = mixed_layer
self.box = box
self.min_level = min_level
self.max_level = max_level
def __eq__(self, other):
if self._different_type(other):
return False
return self.startdate == other.startdate and self.member == other.member and self.chunk == other.chunk and \
self.box == other.box and self.basin == other.basin and self.mxloption == other.mxloption
def __str__(self):
return 'Heat content Startdate: {0.startdate} Member: {0.member} Chunk: {0.chunk} Mixed layer: {0.mxloption} ' \
'Box: {0.box} Basin: {0.basin}'.format(self)
[docs] @classmethod
def generate_jobs(cls, diags, options):
"""
Create a job for each chunk to compute the diagnostic
:param diags: Diagnostics manager class
:type diags: Diags
:param options: basin, mixed layer option (1 to only compute at the mixed layer, -1 to exclude it, 0 to ignore),
minimum depth, maximum depth
:type options: list[str]
:return:
"""
options_available = (DiagnosticBasinOption('basin'),
DiagnosticIntOption('mixed_layer', None, -1, 1),
DiagnosticIntOption('min_depth'),
DiagnosticIntOption('max_depth'))
options = cls.process_options(options, options_available)
box = Box()
box.min_depth = options['min_depth']
box.max_depth = options['max_depth']
min_level = 0
max_level = 0
if box.min_depth > 0 or box.max_depth > 0:
max_level, min_level = cls._get_levels_from_meters(box)
job_list = list()
for startdate, member, chunk in diags.config.experiment.get_chunk_list():
job_list.append(HeatContent(diags.data_manager, startdate, member, chunk,
options['basin'], options['mixed_layer'], box, min_level, max_level))
return job_list
@classmethod
def _get_levels_from_meters(cls, box):
depth_t, depth_w = cls._read_level_values()
def find_nearest(array, value):
idx = (np.abs(array - value)).argmin()
return idx
if box.min_depth > 0:
min_level = find_nearest(depth_t, box.min_depth)
if depth_t[min_level] < box.min_depth:
min_level += 1
else:
min_level = 0
numlevels = len(depth_t)
if box.max_depth > 0:
max_level = find_nearest(depth_t, box.max_depth)
if depth_t[max_level] >= box.max_depth:
max_level -= 1
else:
max_level = numlevels - 1
if min_level < 0:
min_level = 0
if max_level >= numlevels:
max_level = numlevels - 1
box.min_depth = round(depth_w[min_level])
if max_level < numlevels - 1:
box.max_depth = round(depth_w[max_level + 1])
else:
box.max_depth = round(depth_w[-1])
min_level += 1
max_level += 1
return max_level, min_level
@classmethod
def _read_level_values(cls):
handler = Utils.open_cdf('mesh_zgr.nc')
if 'gdepw_1d' in handler.variables:
depth_w = handler.variables['gdepw_1d'][0, :]
elif 'gdepw_0' in handler.variables:
depth_w = handler.variables['gdepw_0'][0, :]
else:
raise Exception('gdepw 1D variable can not be found')
if 'gdept_1d' in handler.variables:
depth_t = handler.variables['gdept_1d'][0, :]
elif 'gdept_0' in handler.variables:
depth_t = handler.variables['gdept_0'][0, :]
else:
raise Exception('gdept 1D variable can not be found')
handler.close()
return depth_t, depth_w
[docs] def request_data(self):
"""Request data required by the diagnostic"""
self.thetao = self.request_chunk(ModelingRealms.ocean, 'thetao', self.startdate, self.member, self.chunk)
if self.mxloption != 0:
self.mlotst = self.request_chunk(ModelingRealms.ocean, 'mlotst', self.startdate, self.member, self.chunk)
[docs] def declare_data_generated(self):
"""Declare data to be generated by the diagnostic"""
if self.box.min_depth == 0:
# For cdftools, this is all levels
box_save = None
else:
box_save = self.box
self.heatcsum = self.declare_chunk(ModelingRealms.ocean, 'heatcsum', self.startdate, self.member, self.chunk,
box=box_save, region=self.basin.name)
self.heatcmean = self.declare_chunk(ModelingRealms.ocean, 'heatcvmean', self.startdate, self.member, self.chunk,
box=box_save, region=self.basin.name)
[docs] def compute(self):
"""Run the diagnostic"""
nco = Utils.nco
temperature_file = TempFile.get()
Utils.copy_file(self.thetao.local_file, temperature_file)
if self.mxloption != 0:
nco.ncks(input=self.mlotst.local_file, output=temperature_file, options=('-A -v mlotst',))
para = list()
if self.min_level != 0:
para.append('-zoom')
para.append(0)
para.append(0)
para.append(0)
para.append(0)
para.append(self.min_level)
para.append(self.max_level)
if self.mxloption != 0:
para.append('-mxloption')
para.append(str(self.mxloption))
if self.basin != Basins().Global:
handler = Utils.open_cdf('mask_regions.3d.nc')
if self.basin.name not in handler.variables:
raise Exception('Basin {0} is not defined on mask_regions.nc'.format(self.basin.name))
handler.close()
para.append('-M')
para.append('mask_regions.3d.nc')
para.append(self.basin.name)
temp2 = TempFile.get()
cdftools.run('cdfheatc', options=para, input_file=temperature_file, output_file=temp2, input_option='-f')
results = Utils.open_cdf(temp2)
heatcsum_temp = TempFile.get()
heatcvmean_temp = TempFile.get()
nco.ncks(input=temperature_file, output=heatcsum_temp, options=('-O -v time',))
shutil.copy(heatcsum_temp, heatcvmean_temp)
heatcsum_handler = Utils.open_cdf(heatcsum_temp)
thc = heatcsum_handler.createVariable('heatcsum', float, 'time')
thc.standard_name = "integral_of_sea_water_potential_temperature_expressed_as_heat_content"
thc.long_name = "Total heat content"
thc.units = "J"
thc[:] = results.variables['heatc3d'][:, 0, 0]
heatcsum_handler.close()
heatcvmean_handler = Utils.open_cdf(heatcvmean_temp)
uhc = heatcvmean_handler.createVariable('heatcvmean', float, 'time')
uhc.standard_name = "integral_of_sea_water_potential_temperature_expressed_as_heat_content"
uhc.long_name = "Heat content per unit volume"
uhc.units = "J*m^-3"
uhc[:] = results.variables['heatc3dpervol'][:, 0, 0]
heatcvmean_handler.close()
results.close()
Utils.setminmax(heatcsum_temp, 'heatcsum')
self.heatcsum.set_local_file(heatcsum_temp)
Utils.setminmax(heatcvmean_temp, 'heatcvmean')
self.heatcmean.set_local_file(heatcvmean_temp)