# coding=utf-8 import numpy as np from box import Box from diagnostic import Diagnostic from earthdiagnostics import cdo from utils import Utils, TempFile class HeatContentLayer(Diagnostic): """ Point-wise Ocean Heat Content in a specified ocean thickness (J/m-2) :original author: Isabel Andreu Burillo :contributor: Virginie Guemas :contributor: Javier Vegas-Regidor :created: June 2012 :last modified: June 2016 """ def __init__(self, data_manager, startdate, member, chunk, box): Diagnostic.__init__(self, data_manager) self.startdate = startdate self.member = member self.chunk = chunk self.box = box self.required_vars = ['so', 'mlotst'] self.generated_vars = ['scvertsum'] def __str__(self): return 'Heat content layer Startdate: {0} Member: {1} Chunk: {2}'.format(self.startdate, self.member, self.chunk) @classmethod def generate_jobs(cls, diags, options): """ Creates a job for each chunk to compute the diagnostic :param diags: Diagnostics manager class :type diags: Diags :param options: None :type options: list[str] :return: """ num_options = len(options) - 1 if num_options < 2: raise Exception('You must specify the minimum and maximum depth to use') if num_options > 2: raise Exception('You must specify 2 for the heat content layer diagnostic') box = Box(True) box.min_depth = int(options[1]) box.max_depth = int(options[2]) job_list = list() for startdate, member, chunk in diags.get_chunk_list(): job_list.append(HeatContentLayer(diags.data_manager, startdate, member, chunk, box)) return job_list def compute(self): nco = Utils.nco temp = TempFile.get() heatc_sl_out = TempFile.get() heatc_sl_top = TempFile.get() level_above = TempFile.get() level_below = TempFile.get() heatc_sl_bottom = TempFile.get() heatc_sl_top_invert = TempFile.get() e3tfile = TempFile.get() results = TempFile.get() handler = Utils.openCdf('mesh_zgr.nc') if 'e3t' in handler.variables: e3t_name = 'e3t' elif 'e3t_0' in handler.variables: e3t_name = 'e3t_0' else: raise Exception('e3t variable can not be found') handler.close() nco.ncap2(input='mesh_zgr.nc', output=e3tfile, options='-v -O -s "heatc_sl=tmask*{0}"'.format(e3t_name)), thetao_file = self.data_manager.get_file('ocean', 'thetao', self.startdate, self.member, self.chunk) nco.ncks(input=thetao_file, output=temp, options='-O -v thetao') Utils.rename_variable(temp, 'thetao', 'heatc_sl') cdo.mul(input=' '.join([temp, e3tfile]), output=heatc_sl_out) # extract the data between the two given depths --> heatc_sl_top.nc nco.ncks(input=heatc_sl_out, output=heatc_sl_top, options='-O -d lev,{0}.0,{1}.0'.format(self.box.min_depth, self.box.max_depth)) # now extract a few levels below, to compute the residual ohc # addition with float returned: nco.ncks(input=heatc_sl_out, output=heatc_sl_bottom, options='-O -d lev,{0}.0,{1}.0'.format(self.box.max_depth, self.box.max_depth + 200)) # obtain the weight for the extra level containing the 300 m # deptht in the gridT files is positive # weight = (300.0 - depth_top)/(depth_bottom - depth_top) # and add the thickness down to 300 m in the next layer nco.ncpdq(options="-a '-lev'", input=heatc_sl_top, output=heatc_sl_top_invert) nco.ncks(input=heatc_sl_top_invert, output=level_above, options='-O -d lev,0,0,1') nco.ncks(input=heatc_sl_bottom, output=level_below, options='-O -d lev,0,0,1') handler = Utils.openCdf(level_above) lev_above = handler.variables['lev'][:] handler.close() handler = Utils.openCdf(level_below) layerthcknss = handler.variables['lev'][:] - lev_above heatc_sl_below = handler.variables['heatc_sl'][:] handler.close() factor = (self.box.max_depth - lev_above) / layerthcknss heatc_sl_below = heatc_sl_below * factor handler = Utils.openCdf(heatc_sl_top) heatc_sl = handler.variables['heatc_sl'][:] handler.close() heatc_sl = np.sum(heatc_sl, 1) heatc_sl = heatc_sl[:] + heatc_sl_below[:, 0, :] heatc_sl = heatc_sl[:] * 1020 * 4000 nco.ncks(input=thetao_file, output=results, options='-O -v lon,lat,time') Utils.rename_variables(results, {'x': 'i', 'y': 'j'}, False, True) handler_results = Utils.openCdf(results) handler_results.createVariable('ohc', float, ('time', 'j', 'i')) handler_results.close() handler_results = Utils.openCdf(results) handler_results.variables['ohc'][:] = heatc_sl handler_results.close() Utils.setminmax(results, 'ohc') self.data_manager.send_file(results, 'ocean', 'ohc', self.startdate, self.member, self.chunk, box=self.box)