# Distribution correction
## A python script to correct distributions using Fortran for speeding up

# Fortran compilation of three modules: module_definintions.f90,  module_generic.f90, 
#   module_scientific.f90 and module_DistricCorrection.f90
# $ make clean
# $ make

# L. Fita, CIMA July 2017

## DistriCorrection.py -f [file] -m [method] -S [values] -d [DistriDimns] -v [varname] -D [debug]
# [file]: file to use
# [method]: method to use
# [values]: specific values for the method
# [DistriDimns]: ',' list of dimensions of the variable to take to create the Distribution
# [varname]: name of the variable to use
# [debug]: debugging flag (False, default)

# Methods _______
## 'randLowRemovePercent_HighEqui'
### This method, removes values from the lower end of the distribution and 
###   equi-distribute the total amount on a percentile on the high end
###   values = [rangeVals]:[PercentRange]:[KindRemove]:[Nquants]:[Hquant]:[Repval]
###     [rangeVals]: ',' list of values from which establish the ranges of removals
###     [PercentRange]: ',' list of percentages of removing to apply to each range (one value less than [rangeVals]) 
###     [KindRemove]: ',' kind of methodology to select which values to remove
###       'rand': randomly selection within the correspondent range
###     [Nquants]: quantity of percentiles to compute from the distribution
###     [Hquant]: percentile of the distribution of values above which total removed amount will be equi-distributed
###     [Repval]: replacing value to be assigned at the removed values
###       'min': initial minimal value of the distribution
###       'mean': initial mean value of the distribution
###       '[value]': any given value otherwise

## e.g. # DistriCorrection_for.py -f E2OFD_SAM_Rainf_2008_2014_mm3hh.nc -S 0,0.0625,0.09375,0.125:0.8,0.6,0.4:rand:20:0.9:0. -m 'randLowRemovePercent_HighEqui' -d time -v Rainf -D True
### removing 80% of interval of values (0,0.0625]
### removing 60% of interval of values (0.0625,0.09375]
### removing 40% of interval of values (0.09375,0.125]
### increasing the values to the 90 percentile and leaving a 0. on the removed points

from optparse import OptionParser
import numpy as np
from netCDF4 import Dataset as NetCDFFile
import os
import re
import numpy.ma as ma
import subprocess as sub
import random
import module_ForDistriCorrect as fDistCorr

main = 'DistriCorrection_for'

# List of available methods
methods = ['randLowRemovePercent_HighEqui']

# output file name
ofilen = main + '.nc'

# Starting from the beginning
inislice = 0

# Functions (retrieved from PyNCplot, or some news)
##
# add_global_CIMA: Function to add global attributes from 'CIMA' to a given netCDF
# basicvardef: Function to give the basic attributes to a variable
# check_arguments: Function to check the number of arguments if they are coincident
# variable_inf: Class to provide the information from a given variable
# index_vec: Function to provide the coordinates of a given value inside a vector
# get_Variables: Function to get a list of variables from an existing file
# multi_index_mat: Function to provide the multiple coordinates of a given value inside a matrix
# Nstrings: Function to provide repeated N-times a given string
# PrinQuantiles: Function to print the quantiles of values
# provide_slices: Function to provide a list of slices for a matrix giving a 
#   sub-section of running dimensions
# Quantiles: Class to provide quantiles from a given array of values
# searchInlist: Function to search a value within a list
# set_attribute: Sets a value of an attribute of a netCDF variable. Removes previous 
#   attribute value if exists
# set_attributek: Sets a value of an attribute of a netCDF variable with a kind. 
#   Removes previous attribute value if exists
# Str_Bool: Function to transform from a String value to a boolean one
# str_list: Function to obtain a list from a string giving a split character
# str_list_k: Function to obtain a list of types of values from a string giving a split character
# typemod: Function to give back a value in a given dtype
# variable_inf: Class to provide the information from a given variable

errormsg = 'ERROR -- error -- ERROR -- error'
warnmsg = 'WARNING -- warning -- WARNING -- warning'
infomsg = 'INFORMATION -- information -- INFORMATION -- information'

# masked array (for type tests)
mamat = ma.masked_equal([0,1],1)

####### ###### ##### #### ### ## #

def searchInlist(listname, nameFind):
    """ Function to search a value within a list
    listname = list
    nameFind = value to find
    >>> searInlist(['1', '2', '3', '5'], '5')
    True
    """
    for x in listname:
      if x == nameFind:
        return True
    return False

def check_arguments(funcname,args,expectargs,char):
    """ Function to check the number of arguments if they are coincident
    check_arguments(funcname,Nargs,args,char)
      funcname= name of the function/program to check
      args= passed arguments
      expectargs= expected arguments
      char= character used to split the arguments
    """

    fname = 'check_arguments'

    Nvals = len(args.split(char))
    Nargs = len(expectargs.split(char))

    if Nvals != Nargs:
        print errormsg
        print '  ' + fname + ': wrong number of arguments:',Nvals," passed to  '",   \
          funcname, "' which requires:",Nargs,'!!'
        print '     # given expected _______'
        Nmax = np.max([Nvals, Nargs])
        for ia in range(Nmax):
            if ia > Nvals-1:
                aval = ' '
            else:
                aval = args.split(char)[ia]
            if ia > Nargs-1:
                aexp = ' '
            else:
                aexp = expectargs.split(char)[ia]

            print '    ', ia, aval, aexp
        quit(-1)

    return

def index_vec(vec,val):
    """ Function to provide the coordinates of a given value inside a vector
    index_vec(vec,val)
      vec= vector with values
      val= value to search
    >>> index_vec(np.arange(27),22)
    22
    """

    fname = 'index_vec'

    vecv = np.array(vec)

    valpos = -1
    for i in range(vecv.shape[0]):
        if vecv[i] == val:
            valpos = i
            break

    return valpos

def provide_slices(dimns, dimzs, rundims):
    """ Function to provide a list of slices for a matrix giving a sub-section of running dimensions
        dimns: list of names of dimensions
        dimzs: list of sizes of dimensions in the same order as in [dimns]
        rundims: name of dimensions to allow to run
    >>> provide_slices(['z', 't', 'y', 'x'], [2, 3, 10, 20], ['t', 'z', 'l'])
    [[1, 2, slice(0, 10, None), slice(0, 20, None)], 
     [0, 1, slice(0, 10, None), slice(0, 20, None)], 
     [1, 1, slice(0, 10, None), slice(0, 20, None)], 
     [0, 0, slice(0, 10, None), slice(0, 20, None)], 
     [1, 0, slice(0, 10, None), slice(0, 20, None)], 
     [0, 2, slice(0, 10, None), slice(0, 20, None)]]
    """
    fname = 'provide_slices'

    Ndims = len(dimns)

    # Checking presence of dimension in list of dimensions
    Tslicesize = 1
    origslicesize = {}
    runslicesize = {}
    for dimn in rundims:
        if not searchInlist(dimns, dimn):
            #print warnmsg
            #print '  ' + fname + ": dimension '" + dimn + "' not in list !!"
            #print '    removing it'
            rundims.remove(dimn)

    for dimn in dimns:
        idim = index_vec(dimns, dimn)
        origslicesize[dimn] = dimzs[idim]
        if searchInlist(rundims, dimn):
            runslicesize[dimn] = dimzs[idim] - 1
            Tslicesize = Tslicesize*dimzs[idim]

    #print '  ' + fname + ': Total number of slices to provide:', Tslicesize,         \
    #  'along:', runslicesize

    slices = []
    for il in range(Tslicesize):
        islice = []
        alreadychanged = False
        for idim in range(Ndims):
            dimn = dimns[idim]
            if searchInlist(rundims,dimn):
                islice.append(runslicesize[dimn])
                # Only running value for a given dimension for the next
                if not alreadychanged:
                    runslicesize[dimn] = runslicesize[dimn] - 1
                    alreadychanged = True
                    if runslicesize[dimn] < 0:
                        runslicesize[dimn] = origslicesize[dimn] - 1
                        for rdimn in rundims:
                            if rdimn != dimn:
                                runslicesize[rdimn] = runslicesize[rdimn] - 1
                                if runslicesize[rdimn] >= 0:
                                    break
                                else:
                                    runslicesize[rdimn] = origslicesize[rdimn] - 1
            else:
                islice.append(slice(0,origslicesize[dimn]))

        slices.append(islice)

    return slices

def str_list(string, cdiv):
    """ Function to obtain a list from a string givin a split character
      string= String from which to obtain a list ('None' for None)
      cdiv= character to use to split the string
    >>> str_list('1:@#:$:56', ':')
    ['1', '@#', '$', '56']
    >>> str_list('1@#$56', ':')
    ['1@#$56']
    """
    fname = 'str_list'

    if string.find(cdiv) != -1:
        listv = string.split(cdiv)
    else:
        if string == 'None':
            listv = None
        else:
            listv = [string]

    return listv

def typemod(value, typeval):
    """ Function to give back a value in a given dtype
    >>> print(typemod(8.2223, 'np.float64'))
    <type 'numpy.float64'>
    >>> print(typemod(8.2223, 'tuple'))
    <type 'tuple'>
    """
    fname='typemod'

    if typeval == 'int':
        return int(value)
    elif typeval == 'long':
        return long(value)
    elif typeval == 'float':
        return float(value)
    elif typeval == 'complex':
        return complex(value)
    elif typeval == 'str' or typeval == 'S':
        return str(value)
    elif typeval == 'bool':
        return bool(value)
    elif typeval == 'list':
        newval = []
        newval.append(value)
        return newval
    elif typeval == 'dic':
        newval = {}
        newval[value] = value
        return newval
    elif typeval == 'tuple':
        newv = []
        newv.append(value)
        newval = tuple(newv)
        return newval
    elif typeval == 'np.int8':
        return np.int8(value)
    elif typeval == 'np.int16':
        return np.int16(value)
    elif typeval == 'np.int32':
        return np.int32(value)
    elif typeval == 'np.int64':
        return np.int64(value)
    elif typeval == 'np.uint8':
        return np.uint8(value)
    elif typeval == 'np.uint16':
        return np.uint16(value)
    elif typeval == 'np.np.uint32':
        return np.uint32(value)
    elif typeval == 'np.uint64':
        return np.uint64(value)
    elif typeval == 'np.float' or typeval == 'R':
        return np.float(value)
    elif typeval == 'np.float16':
        return np.float16(value)
    elif typeval == 'np.float32':
        return np.float32(value)
    elif typeval == 'float32':
        return np.float32(value)
    elif typeval == 'np.float64':
        return np.float64(value)
    elif typeval == 'float64':
        return np.float64(value)
    elif typeval == 'np.complex':
        return np.complex(value)
    elif typeval == 'np.complex64':
        return np.complex64(value)
    elif typeval == 'np.complex128':
        return np.complex128(value)
    else:
        print errormsg
        print fname + ':  data type "' + typeval + '" is not ready !'
        print errormsg
        quit(-1)

    return

def str_list_k(string, cdiv, kind):
    """ Function to obtain a list of types of values from a string giving a split character
      string= String from which to obtain a list ('None' for None)
      cdiv= character to use to split the string
      kind= kind of desired value (as string like: 'np.float', 'int', 'np.float64', ....)
    >>> str_list_k('1:@#:$:56', ':', 'S')
    ['1', '@#', '$', '56']
    >>> str_list_k('1:3.4:12.3', ':', 'np.float64')
    [1.0, 3.3999999999999999, 12.300000000000001]
    """
    fname = 'str_list'

    if string.find(cdiv) != -1:
        listv = string.split(cdiv)
    else:
        if string == 'None':
            listv = None
        else:
            listv = [string]

    if listv is not None:
        finalist = []
        for lv in listv:
            finalist.append(typemod(lv, kind))
    else:
        finalist = None

    return finalist

class Quantiles(object):
    """ Class to provide quantiles from a given array of values
    """

    def __init__(self, values, Nquants):
        import numpy.ma as ma

        if values is None:
            self.quantilesv = None

        else:
            self.quantilesv = []

            vals0 = values.flatten()
            Nvalues = len(vals0)
            vals = ma.masked_equal(vals0, None)
            Nvals=len(vals.compressed())

            sortedvals = sorted(vals.compressed())
            for iq in range(Nquants):
                self.quantilesv.append(sortedvals[int((Nvals-1)*iq/Nquants)])

            self.quantilesv.append(sortedvals[Nvals-1])

        return

def Str_Bool(val):
    """ Function to transform from a String value to a boolean one
    >>> Str_Bool('True')
    True
    >>> Str_Bool('0')
    False
    >>> Str_Bool('no')
    False
    """

    fname = 'Str_Bool'

    truelist = ['True', 'true', '1', 'Yes', 'yes', 'T', 't'] 
    falselist = ['False', 'false', '0', 'No', 'no', 'F', 'f']

    if searchInlist(truelist,val): 
        boolv = True
    elif searchInlist(falselist, val):
        boolv = False
    else:
        print errormsg
        print '  ' + fname + ": value '" + val + "' not ready!!"
        quit(-1)

    return boolv

def multi_index_mat(mat,val):
    """ Function to provide the multiple coordinates of a given value inside a matrix
    index_mat(mat,val)
      mat= matrix with values
      val= value to search
    >>> vals = np.ones((24), dtype=np.float).reshape(2,3,4)
    vals[:,:,2] = 0.
    vals[1,:,:] = np.pi
    vals[:,2,:] = -1.
    multi_index_mat(vals,1.)
    [array([0, 0, 0]), array([0, 0, 1]), array([0, 0, 3]), array([0, 1, 0]), array([0, 1, 1]), array([0, 1, 3])]
    """
    fname = 'multi_index_mat'

    matshape = mat.shape

    ivalpos = []
    matlist = list(mat.flatten())
    Lmatlist = len(matlist)
    
    val0 = val - val
    if val != val0:
        valdiff = val0
    else:
        valdiff = np.ones((1), dtype = type(val))
    
    ifound = 0
    while ifound < Lmatlist:
        if matlist.count(val) == 0:
            ifound = Lmatlist + 1
        else:
            ifound = matlist.index(val)

            Ndims = len(matshape)
            valpos = np.zeros((Ndims), dtype=int)
            baseprevdims = np.zeros((Ndims), dtype=int)

            for dimid in range(Ndims):
                baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims])
                if dimid == 0:
                    alreadyplaced = 0
                else:
                    alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
                valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
            matlist[ifound] = valdiff
            ivalpos.append(valpos)

    return ivalpos

def Nstrings(strv, Ntimes):
    """ Function to provide repeated N-times a given string
      strv: string value
      Ntimes: times to repeat the string
    >>> Nstrings('123', 4)
    123123123123
    """
    fname = 'Nstrings'

    newstring = ''
    for i in range(Ntimes):
        newstring = newstring + strv

    return newstring

def PrinQuantiles(vals, qtvs, Nspcs=6):
    """ Function to print the quantiles of values
      vals: series of values
      qtvs: values of the quantiles
      Nspcs: base quantity of spaces
    """
    fname = 'PrinQuantiles'

    # There is Nqts + 1
    Nqts = len(qtvs)

    bspc = Nstrings(' ',Nspcs)
    for iq in range(Nqts-1):
        search1 = vals >= qtvs[iq]
        search2 = vals <= qtvs[iq+1]
        search = search1*search2
        print bspc, '[',iq+1,']','{:6.2f}'.format((iq+1)*100./(Nqts-1)), '%:',       \
          qtvs[iq+1], 'N:', np.sum(search)

    return

####### ###### ##### #### ### ## # netCDF functions # ## ### #### ##### ###### #######

def set_attribute(ncv, attrname, attrvalue):
    """ Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists
    ncv = object netcdf variable
    attrname = name of the attribute
    attrvalue = value of the attribute
    """

    attvar = ncv.ncattrs()
    if searchInlist(attvar, attrname):
        attr = ncv.delncattr(attrname)

    attr = ncv.setncattr(attrname, attrvalue)

    return attr

def set_attributek(ncv, attrname, attrval, attrkind):
    """ Sets a value of an attribute of a netCDF variable with a kind. Removes previous attribute value if exists
    ncvar = object netcdf variable
    attrname = name of the attribute
    attrvalue = value of the attribute
    atrtrkind = kind of attribute: 'S', string ('!' as spaces); 'U', unicode ('!' as spaces); 'I', integer; 
      'Inp32', numpy integer 32; 'R', real; ' npfloat', np.float; 'D', np.float64
    """
    fname = 'set_attributek'
    validk = {'S': 'string', 'U': 'unicode', 'I': 'integer',                         \
      'Inp32': 'integer long (np.int32)', 'R': 'float', 'npfloat': 'np.float',       \
      'D': 'float long (np.float64)'}

    if type(attrkind) == type('s'):
        if attrkind == 'S':
            attrvalue = str(attrval.replace('!', ' '))
        elif attrkind == 'U':
            attrvalue = unicode(attrval.replace('!',' '))
        elif attrkind == 'I':
            attrvalue = np.int(attrval)
        elif attrkind == 'R':
            attrvalue = float(attrval)
        elif attrkind == 'npfloat':
            attrvalue = np.float(attrval)
        elif attrkind == 'D':
            attrvalue = np.float64(attrval)
        else:
            print errormsg
            print '  ' + fname + ": '" + attrkind + "' kind of attribute is not ready!"
            print '  valid values: _______'
            for key in validk.keys():
                print '    ', key,':', validk[key]
            quit(-1)
    else:
        if attrkind == type(str('a')):
            attrvalue = str(attrval.replace('!', ' '))
        elif attrkind == type(unicode('a')):
            attrvalue = unicode(attrval.replace('!',' '))
        elif attrkind == type(np.int(1)):
            attrvalue = np.int(attrval)
        elif attrkind == np.dtype('i'):
            attrvalue = np.int32(attrval)
        elif attrkind == type(float(1.)):
            attrvalue = float(attrval)
        elif attrkind == type(np.float(1.)):
            attrvalue = np.float(attrval)
        elif attrkind == np.dtype('float32'):
            attrvalue = np.array(attrval, dtype='f')
        elif attrkind == type(np.float32(1.)):
            attrvalue = np.float32(attrval)
        elif attrkind == type(np.float64(1.)):
            attrvalue = np.float64(attrval)
        elif attrkind == type(np.array([1.,2.])):
            attrvalue = attrval
        else:
            print errormsg
            print '  ' + fname + ": '" + attrkind + "' kind of attribute is not ready!"
            print '  valid values: _______'
            for key in validk.keys():
                print '    ', key,':', validk[key]
            quit(-1)

    attvar = ncv.ncattrs()
    if searchInlist(attvar, attrname):
        attr = ncv.delncattr(attrname)

    if attrname == 'original_subroutines_author': 
        attrvalue = 'Cindy Bruyere'
    attr = ncv.setncattr(attrname, attrvalue)
        
    return attr

class variable_inf(object):
    """ Class to provide the information from a given variable
    var = object netCDF variable
    self.name: name of the variable
    self.dtype: type of the variable
    self.attributes: list with the name of attributes
    self.FillValue: value of the missing value
    self.dimns: name of the dimensions of the variable
    self.dims: dimensions of the variable
    self.Ndims: number of dimensions
    self.dimx: length of dimension along x-axis
    self.dimy: length of dimension along y-axis
    self.sname: standard name
    self.lname: long name
    self.corr: attribute 'coordinates'
    self.units: units of the variable
    """

    def __init__(self, var):

        if var is None:
            self.name = None
            self.dimns = None
            self.dims = None
            self.Ndims = None
            self.dimx = None
            self.dimy = None
            self.sname = None
            self.lname = None
            self.corr = None
            self.units = None
            self.FillValue = None
            self.dtype = None
            self.attributes = None
        else:
            self.name = var._name
            self.dimns = var.dimensions
            self.dtype = var.dtype
            self.attributes = var.ncattrs()
            self.dims = var.shape

            if gen.searchInlist(self.attributes, 'standard_name'):
                self.sname = var.getncattr('standard_name')
            else:
                self.sname = None

            if gen.searchInlist(self.attributes, 'long_name'):
                self.lname = var.getncattr('long_name')
            else:
                self.lname = None

            if gen.searchInlist(self.attributes, 'coordinates'):
                self.coor = var.getncattr('coordinates')
            else:
                self.coor = None

            if gen.searchInlist(self.attributes, 'units'):
                self.units = var.getncattr('units')
            else:
                self.units = None

            if gen.searchInlist(self.attributes, '_FillValue'):
                self.FillValue = var.getncattr('_FillValue')
            else:
                self.FillValue = None
             
            self.Ndims = len(self.dims)
            if self.Ndims == 1:
                self.dimx=self.dims[0]
            if self.Ndims == 2:
                self.dimy=self.dims[0]
                self.dimx=self.dims[1]
            if self.Ndims == 3:
                self.dimy=self.dims[1]
                self.dimx=self.dims[2]
            if self.Ndims == 4:
                self.dimy=self.dims[2]
                self.dimx=self.dims[3]

        return

def basicvardef(varobj, vstname, vlname, vunits):
    """ Function to give the basic attributes to a variable
    varobj= netCDF variable object
    vstname= standard name of the variable
    vlname= long name of the variable
    vunits= units of the variable
    """
    attr = varobj.setncattr('standard_name', vstname)
    attr = varobj.setncattr('long_name', vlname)
    attr = varobj.setncattr('units', vunits)

    return

def get_Variables(values, ncfile, variables):
    """ Function to get a list of variables from an existing file
      novar= behaviour when no variable is found
        'error': an error messages arise and program stops
        'warn': a warning messages arise and program continues
      ncfile= WRF file to use
      variables = ',' separated list of names of variables to get
    """
    fname = 'get_Variables'

    if values == 'h':
        print fname + '_____________________________________________________________'
        print get_Variables.__doc__
        quit()

    novar = values
    variables = variables.split(',')

    ofile = 'get_Variables.nc'

    onewnc = NetCDFFile(ofile, 'w')

    onc = NetCDFFile(ncfile, 'r')
    filevars = list(onc.variables.keys())

    for var in variables:
        if not gen.searchInlist(filevars, var):
            if novar == 'error':
                print errormsg
                print '  ' + fname + ": file '" + ncfile + "' does not have " +      \
                  "variable '" + var + "' !!"
                print '    variables in file:', filevars
                quit(-1)
            elif novar == 'warn':
                print warnmsg
                print '  ' + fname + ": file '" + ncfile + "' does not have " +      \
                  "variable '" + var + "' !!"
                print '    variables in file:', filevars
                continue
            else:
                print errormsg
                print '  ' + fname + ": novar '" + novar + "' not ready !!"
                quit(-1)

        print '  ' + fname + ": getting variable '" + var + "'..."
        ovar = onc.variables[var]
        varinf = variable_inf(ovar)
        vardims = varinf.dimns

        newfiledims = list(onewnc.dimensions.keys())

        for vdim in vardims:
            if not gen.searchInlist(newfiledims, vdim):
                print '  ' + fname + ": including dimension '" + vdim + "'..."
                if onc.dimensions[vdim].isunlimited():
                    newdim = onewnc.createDimension(vdim, None)
                else:
                    newdim = onewnc.createDimension(vdim, len(onc.dimensions[vdim]))

        if varinf.FillValue is not None:
            newvar = onewnc.createVariable(var, nctype(varinf.dtype), tuple(vardims),\
              fillvalue = varinf.FillValue) 
        else:
            newvar = onewnc.createVariable(var, nctype(varinf.dtype), tuple(vardims)) 
        newvar[:] = ovar[:]
        for attrn in varinf.attributes:
            attrv = ovar.getncattr(attrn)
            newattr = set_attribute(newvar, attrn, attrv)

        onewnc.sync()

    # Global attributes
    for attrn in onc.ncattrs():
        attrv = onc.getncattr(attrn)
        newattr = set_attribute(onewnc, attrn, attrv)
    
    onewnc.sync()
    onc.close()

    onewnc.close()

    print fname + ": successful written of '" + ofile + "' !!"

def add_global_CIMA(ObjFile, pyscript, version):
    """ Function to add global attributes from 'CIMA' to a given netCDF
      ObjFile= netCDF file object to which add the global attributes
      version= version of the function
    """
    fname = 'add_global_PyNCplot'

    # Global values
    ObjFile.setncattr('author', 'L. Fita')
    newattr = set_attributek(ObjFile, 'institution', unicode('Centro de ' +          \
     'Investigaciones del Mar y la Atm' + unichr(243) + 'sfera (CIMA)'), 'U')
    ObjFile.setncattr('university', 'Pierre et Marie Curie - Jussieu, Paris 06')
    ObjFile.setncattr('center', 'Consejo Nacional de Investigaciones Cient' +        \
      unicode(237) + 'ficas y T' + unicode(233) + 'cnicas (CONICET)')
    ObjFile.setncattr('city', 'Buenos Aires')
    ObjFile.setncattr('country', 'Argentina')
    ObjFile.setncattr('tool', 'DistriCorection')
    ObjFile.setncattr('url', 'http://www.xn--llusfb-5va.cat/python/PyNCplot/')
    ObjFile.setncattr('script', pyscript)
    #ObjFile.setncattr('function', funcname)
    ObjFile.setncattr('version', version)

    ObjFile.sync()  

    return

# testing with 100 gaussian values
#exvalues = np.random.normal(12., 4., 100)

####### ###### ##### #### ### ## #

# Arguments
##
parser = OptionParser()
parser.add_option("-d", "--Dimension", dest="dimns", 
  help="',' list of names of the dimension to construct distribution", 
  metavar="LABELS")
parser.add_option("-D", "--Debug", dest="debug", 
  help="debug prints",  metavar="BOOL")
parser.add_option("-f", "--netCDF_file", dest="ncfile", help="file to use", 
  metavar="FILE")
parser.add_option("-m", "--Methodology", type='choice', dest="method", 
  choices=methods, help="available methods to modify distribution", metavar="LABEL")
parser.add_option("-S", "--Values", dest="values", 
  help="values for the given methodology", metavar="VALUES")
parser.add_option("-v", "--Variable", dest="varn", 
  help="name of the variable to use", metavar="LABEL")

(opts, args) = parser.parse_args()

#######    #######
## MAIN
    #######

if opts.debug is None:
    print infomsg
    print '  ' + main + ": no debug value provided!!"
    print '    assuming:', False
    debug = False
else:
    debug = Str_Bool(opts.debug)

if not os.path.isfile(opts.ncfile):
      print errormsg
      print '   ' + main + ": File '" + opts.ncfile + "' does not exist !!"
      quit(-1)    

ncf = NetCDFFile(opts.ncfile,'r')

# Getting variable
if not ncf.variables.has_key(opts.varn):
    print errormsg
    print '   ' +main+ ": File '" + opts.ncfile + "' does not content variable '" +\
      opts.varn + "' !!"
    print '    available ones:', ncf.variables.keys()
    ncf.close()
    quit(-1)

ovar = ncf.variables[opts.varn]
# Getting dimensions
dimns = str_list(opts.dimns, ',')
Nvalues = 1

# List of dimensions along which distributions will be modified
rundimns = list(ovar.dimensions)

for dimn in dimns:
    if not searchInlist(list(ovar.dimensions), dimn):
        print errormsg
        print '   ' +main+ ": Variable '" +opts.varn+ "' does not have dimension'" + \
          dimn + "' !!"
        print '    available ones:', ovar.dimensions
        ncf.close()
        quit(-1)
    rundimns.remove(dimn)
    Nvalues = Nvalues * len(ncf.dimensions[dimn])

# We need to assume that given variable will have a large amount of data (at least 
#   along the given dimension)

varvalues = ovar[:]

Ndistri = np.sum(ovar.shape[1:3])

print infomsg
print '  ' + main + ': there are ', Ndistri, " distributions to modify from '" +     \
  opts.varn + "' along dimensions ", dimns, " with: ", Nvalues, 'values'

# File creation as an initial copy of the original file
##
if inislice == 0:
    sout = sub.call('cp ' + opts.ncfile + ' ' + ofilen, shell=True)
    print main + ": successful creation of '" + ofilen + "' !!"

# re-writing the distribution in the new file
onewnc = NetCDFFile(ofilen, 'a')
onewvar = onewnc.variables[opts.varn]

# Getting units of the variable
vunits = onewvar.units

# Getting variable values
## Loop and run along the values of slices being very inefficient? For that program 
##   in Fortran!

# For checking
#slices = []
#slices.append([slice(0, 20456, None), 125, 50])
#Ndistri = 1

# Modifying distriubution 
if opts.method == 'randLowRemovePercent_HighEqui':

    # Getting configuration
    rangevals = str_list_k(opts.values.split(':')[0], ',', 'np.float')
    percenrange = str_list_k(opts.values.split(':')[1], ',', 'np.float')
    kindremove = opts.values.split(':')[2]
    Nquants = int(opts.values.split(':')[3])
    hquant = np.float64(opts.values.split(':')[4])
    repvalS = opts.values.split(':')[5]

    Nranges = len(rangevals)

    # Getting the arguments for the Fortran subroutine
    dim1 = varvalues.shape[2]
    dim2 = varvalues.shape[1]
    dim3 = varvalues.shape[0]
    fillV = np.float64(varvalues.fill_value)

    tvarvals = varvalues.transpose()
    trangevals = np.array(rangevals, dtype=np.float64).transpose()
    tpercenrange = np.array(percenrange, dtype=np.float64).transpose()

    # Computing new distribution      
    tnewdistri, tinredistributed, tamountranges, taboveVal, thigquantind =           \
      fDistCorr.module_districorrection.randlowremovepercent_highequi3d(             \
      dist=tvarvals, rangevals=trangevals, percentrange=tpercenrange,                \
      kindremove=kindremove, nquants=Nquants, hquant=hquant, repvals=repvalS,        \
      fillvalue=fillV, dbg=debug)

    newdistri = tnewdistri.transpose()
    inredistributed = tinredistributed.transpose()
    amountranges = tamountranges.transpose()
    aboveVal = taboveVal.transpose()
    higquantind = thigquantind.transpose()

    onewvar[:] = newdistri
    onewnc.sync()

    # Adding new variables to the output file
    rangevals = str_list_k(opts.values.split(':')[0], ',', 'np.float')
    percenrange = str_list_k(opts.values.split(':')[1], ',', 'np.float')
    kindremove = opts.values.split(':')[2]
    Nquants = int(opts.values.split(':')[3])
    hquant = np.float(opts.values.split(':')[4])
    repvalS = opts.values.split(':')[5]

    Nranges = len(rangevals)

    # Dimensions
    newdim = onewnc.createDimension('range',Nranges-1)
    newdim = onewnc.createDimension('rng_bounds',2)

    # Variable
    #newvar = onewnc.createVariable('methodology','c')
    #newvar[:] = opts.method

    newvar = onewnc.createVariable('ranges','f',('rng_bounds','range'))
    for ir in range(Nranges-1):
        newvar[0,ir] = rangevals[ir]
        newvar[1,ir] = rangevals[ir+1]
    basicvardef(newvar, 'range_values', 'range of values (,] to modify distribution',\
      ovar.units)
    newattr = set_attribute(newvar, 'removing_methodology', kindremove)
    
    newvar = onewnc.createVariable('percentages','f',('range'))
    newvar[:] = percenrange[:]
    basicvardef(newvar,'range_percentages','percentages of removal of each range','1')

    newvar = onewnc.createVariable('amount','f',('range','lat','lon'))
    newvar[:] = amountranges[tuple([slice(0,Nranges-1),slice(0,dim2),slice(0,dim1)])]
    basicvardef(newvar,'range_amounts','amounts of removal of each range', vunits)
    
    newvar = onewnc.createVariable('highquant','f', ('lat', 'lon'))
    newvar[:] = aboveVal
    basicvardef(newvar, 'highquant', 'High quantile above and equal which ' +         \
      'increase values', vunits)

    onewnc.sync()

# Global attributes
add_global_CIMA(onewnc, main, '0.1')
newattr = set_attribute(onewnc, 'modification_method', opts.method)

onewnc.sync()
onewnc.close()