heatcontentlayer.py

# coding=utf-8
import numpy as np

from earthdiagnostics.box import Box
from earthdiagnostics.diagnostic import Diagnostic
from earthdiagnostics import cdo
from earthdiagnostics.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 <virginie.guemas@bsc.es>
    :contributor: Javier Vegas-Regidor<javier.vegas@bsc.es>

    :created: June 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 box: box to use for the calculations
    :type box: Box
    """

    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.exp_manager.get_chunk_list():
            job_list.append(HeatContentLayer(diags.data_manager, startdate, member, chunk, box))
        return job_list

    def compute(self):
        """
        Runs the diagnostic
        """
        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)