# Python script to generate the files and the plots for a sensitivity study with multiple modle configurations and models
#   L. Fita
#   LMD, Jussieu/Palaisseau, France
#  Script configuration get from ASCII file 'model_graphics.dat'

# `Wheeler-Kiladis' diagram: https://github.com/UV-CDAT/wk
#  $ cd /home/lluis/Downloads
#  $ git clone https://github.com/UV-CDAT/wk
#  $ cd wk
#  $ su
#  # python setup.py build >& run_pysetup.log
#  # python setup.py install >& run_pyinstall.log

import numpy as np
import nc_var_tools as ncvar
import generic_tools as gen
import drawing_tools as drw
import time as tim
#  To avoid errors like: '/bin/sh: 1: Syntax error: Bad fd number'
#    Make sure that /bin/sh directs to bash and not to dash: 
#  http://stackoverflow.com/questions/15809060/sh-syntax-error-bad-fd-number
import subprocess as sub
import os

main = 'model_graphics.py'

errmsg = ncvar.errormsg
warnmsg = ncvar.warnmsg

def verify_configuration(ExpConf, debug):
    """ Function to verify mandatory keys from experiment configuration
      ExpConf= dictionary with the configuration
    """
    fname = 'verify_configuration'
 
    # List with the mandatory keys
    mandatorykeys = ['pyHOME', 'cdoHOME', 'scratch', 'filescratch', 'figscratch',    \
      'diffscratch', 'figdiffscratch', 'figallmodexpscratch', 'addfiles',            \
      'addfigures', 'adddiffs', 'adddifffigures', 'addallmodexpfigures', 'debug',    \
      'models', 'modgraphchar', 'expgraphchar', 'ifold', 'ofold', 'warnmsg',         \
      'errmsg', 'titleoperations', 'kindfig', 'CFreftime', 'CFunitstime', 'mapval',  \
      'timekind', 'timefmt', 'timelabel' ]

    # List with the optional keys
    optionalkeys = ['specificvarplot', 'specificdiffopplot', 'specificdiffvarplot',  \
      'RefProj']

    # Dictionary with the mandatory keys which are required as function of a given 
    #   key from the dictionary
    #     if [key] must exist [keyB] = [values]
    ifkeys = {'WRFexps': ['models','WRF'], 'LMDZexps': ['models','LMDZ'],            \
      'WRF_LMDZexps': ['models','WRF_LMDZ'], 'DYNAMICOexps': ['models','DYNAMICO'],  \
      'mthDYNAMICOexps': ['models','mthDYNAMICO'],                                   \
      'WRFheaders': ['models','WRF'], 'WRF_LMDZheaders': ['models','WRF_LMDZ'],      \
      'LMDZheaders': ['models','LMDZ'], 'DYNAMICOheaders': ['models','DYNAMICO'],    \
      'mthDYNAMICOheaders': ['models','mthDYNAMICO']}

    # Dictionary with the optional keys (at least one) which are required as function
    #   of a given key from the dictionary
    #     if [key] must exist at least one of the [keyA, keyB, keyC, ...] with values
    ifopkeys = {'RefProj': ['reprojectvar_remapbil', 'reprojectvar_remapbic',        \
      'reprojectvar_remapdis', 'reprojectvar_remapnn', 'reprojectvar_remapcon',      \
      'reprojectvar_remapcon2', 'reprojectvar_remaplaf', 'reprojectvar_dis',         \
      'reprojecvar_dis']}

    for kdict in mandatorykeys:
        if not ExpConf.has_key(kdict):
            print errmsg
            print '  ' + fname + "' configuration without required '" + kdict + "' !!"
            quit(-1)
        if debug:
            print ExpConf[kdict]

    for kdict in ifkeys.keys():
        vals = ifkeys[kdict]
        keyn = vals[0]
        keyv = vals[1]
        if type(ExpConf[keyn]) == type('A'):
            if not ExpConf.has_key(kdict) and ExpConf[keyn] == keyv:
                print errmsg
                print '  ' + fname + "' configuration without '" + kdict +           \
                  "' when it is required with configuration '" + keyn + '=' + keyv + \
                  "' !!"
                quit(-1)
        elif type(ExpConf[keyn]) == type(['A', 'B']):
            keyvals = ExpConf[keyn]
            if not ExpConf.has_key(kdict) and gen.searchInlist(keyvals, keyv):
                print errmsg
                print '  ' + fname + "' configuration without '" + kdict +           \
                  "' when it is required with configuration '" + keyn + '= [..., ' + \
                 keyv + ", ...]' !!"
                quit(-1)
        else:
            print errmsg
            print '  ' +fname+ 'dictionary type ', type(ExpConf[keyn]), ' not ready!!'
            quit(-1)
        if debug:
            print '  Experiments configuration _______'
            gen.printing_dictionary(ExpConf)

    for kdict in ifopkeys.keys():
        if ExpConf.has_key(kdict):
            opkeys = ifopkeys[kdict]
            hasopkey = False
            opdictvals = {}
            for opkv in opkeys:
                if ExpConf.has_key(opkv):
                    hasopkey = True
                    opdictvals[opkv] = ExpConf[opkv]

            if not hasopkey:
                print errmsg
                print '  ' + fname + "' configuration wit '" + kdict + "' but " +    \
                  "without any key value:", opkeys, "' when it is required !!"
                quit(-1)

            if debug:
                print "  Optional key '" + kdict + "' with values ________"
                gen.printing_dictionary(opdictvals)

    print fname + ': configuration seems to be fine'

    return

def scratches(config):
    """ Function to set-up if it is needed to start from the scratch
      config: dictionary with the configuration
    """
    fname = 'scratches'

    if gen.Str_Bool(config['scratch']):
        scr = True
        print warnmsg
        print "  " + fname + ": starting from the SCRATCH !!"
        print "    10 seconds left!!"
        filescr = True
        figscr = True
        difscr = True
        figdifscr = True
        moddifscr = True
        figmoddifscr = True
        figallmodexpscr = True
        tim.sleep(10)
    else:
        scr = False
        if gen.Str_Bool(config['filescratch']):
            filescr = True
            print warnmsg
            print "  " + main + ": files starting from the SCRATCH !!"
            print "    5 seconds left!!"
            tim.sleep(5)
        else:
            filescr = False

        if gen.Str_Bool(config['figscratch']):
            figscr = True
            print warnmsg
            print "  " + main + ": figures starting from the SCRATCH !!"
            print "    5 seconds left!!"
            tim.sleep(5)
        else:
            figscr = False

        if gen.Str_Bool(config['diffscratch']):
            difscr = True
            print warnmsg
            print "  " + main + ": experiment differences starting from the SCRATCH !!"
            print "    5 seconds left!!"
            tim.sleep(5)
        else:
            difscr = False

        if gen.Str_Bool(config['figdiffscratch']):
            figdifscr = True
            print warnmsg
            print "  " + main + ": figures experiment differences starting from the SCRATCH !!"
            print "    5 seconds left!!"
            tim.sleep(5)
        else:
            figdifscr = False

        if gen.Str_Bool(config['moddiffscratch']):
            moddifscr = True
            print warnmsg
            print "  " + main + ": model differences starting from the SCRATCH !!"
            print "    5 seconds left!!"
            tim.sleep(5)
        else:
            moddifscr = False

        if gen.Str_Bool(config['figmoddiffscratch']):
            figmoddifscr = True
            print warnmsg
            print "  " + main + ": figures model differences starting from the SCRATCH !!"
            print "    5 seconds left!!"
            tim.sleep(5)
        else:
            figmoddifscr = False

        if gen.Str_Bool(config['figallmodexpscratch']):
            figallmodexpscr = True
            print warnmsg
            print "  " + main + ": figures all model-experiment starting from the SCRATCH !!"
            print "    5 seconds left!!"
            tim.sleep(5)
        else:
            figallmodexpscr = False

    if gen.Str_Bool(config['addall']):
        print warnmsg
        print "  " + fname + ": adding new values everywhere !!"
        addfils = True
        addfigs = True
        adddiffs = True
        adddifffigs = True
        addmoddiffs = True
        addmoddifffigs = True
        addallmodexpfigs = True
    else:
        addfils = gen.Str_Bool(config['addfiles'])
        addfigs = gen.Str_Bool(config['addfigures'])
        adddiffs = gen.Str_Bool(config['adddiffs'])
        adddifffigs = gen.Str_Bool(config['adddifffigures'])
        addmoddiffs = gen.Str_Bool(config['addmoddiffs'])
        addmoddifffigs =  gen.Str_Bool(config['addmoddifffigures'])
        addallmodexpfigs = gen.Str_Bool(config['addallmodexpfigures'])

    debug = gen.Str_Bool(config['debug'])

    return scr, filescr, figscr, difscr, figdifscr, moddifscr, figmoddifscr,         \
      figallmodexpscr, addfils, addfigs, adddiffs, adddifffigs, addmoddiffs,         \
      addmoddifffigs, addallmodexpfigs, debug


class gc(object):
    """ Class which holds all the graphical characteristics
      [label]: label as it should appear in the graphics
      [color]: specific color of the model
      [linetype]: specific type of line of the model
      [marker]: specific marker of the model
      [sizes]: line width and point size for the model
      [tmod]: time modification to apply to the model files (None for nothing)
        'setorigin',[YYYYMMDDHHMISS]: re-set origin of times at [YYYYMMDDHHMISS]
    """
    def __init__( self, label, color, linetype, marker, sizes, Tmod):
        self.label = None
        self.color = None
        self.linetype = None
        self.marker = None
        self.sizes = None
        self.tmod = None
        if label is not None:
            self.label = label
            self.color = color
            self.linetype = linetype
            self.marker = marker
            self.sizes = sizes
            if Tmod != 'None':
                self.tmod = Tmod

def graphical_characteristics(config, debug):
    """ Function to provide the graphical characteristics of models and experiments
      config= configuration of the experiment
    """
    fname = 'graphical_characteristics'

    modgraphchar = {}
    modvals = config['modgraphchar'].split(':')
    for modval in modvals:
        modv = modval.split('|')
        modgraphchar[modv[0]] = gc(modv[1], modv[2], modv[3], modv[4], modv[5], modv[6])

    expgraphchar = {}
    expvals = config['expgraphchar'].split(':')
    for expval in expvals:
        expvs = expval.split('|')
        modexp = expvs[0] + '/' + expvs[1]
        if not modgraphchar.has_key(expvs[0]):
            print errmsg
            print '  ' + fname + ": a model called '" +expvs[0]+ "' does not exist !!"
            print '    existing models with graphic characteristics:',               \
              modgraphchar.keys()
            quit(-1)

        modvs = modgraphchar[expvs[0]]
        expv = expvs[2:7]
        if expv[1] == 'asmodel': expv[1] = modvs.color
        if expv[2] == 'asmodel': expv[2] = modvs.linetype
        if expv[3] == 'asmodel': expv[3] = modvs.marker
        if expv[4] == 'asmodel': expv[4] = modvs.sizes

        expgraphchar[modexp] = gc(expv[0], expv[1], expv[2], expv[3], expv[4], 'None')

    # Running some tests
    mods = config['models'].split(':')
    for mod in mods:
        if not modgraphchar.has_key(mod):
            print errmsg
            print '  ' + fname + ": model called '" + mod + "' does not have " +     \
              "graphical characteristics !!"
            print '    existing models with graphic characteristics:',               \
              modgraphchar.keys()
            quit(-1)

        exps = config[mod + 'exps'].split(':')
        for exp in exps:
            modexp = mod + '/' + exp
            if not expgraphchar.has_key(modexp):
                print errmsg
                print '  ' + fname + ": model/experiment called '" + modexp +        \
                  "' does not have graphical characteristics !!"
                print '    existing model/experiments with graphic ' +               \
                  'characteristics:', expgraphchar.keys()
                quit(-1)

    if debug:
        print '  ' + fname + ': model graphical characteristics _______'
        for modn in modgraphchar.keys():
            with gen.Capturing() as output:
                gen.printing_class(modgraphchar[modn])
            print "'" + modn + "'; " + ', '.join(output)
        print '  ' + fname + ': experiment graphical characteristics _______'
        for modexpn in expgraphchar.keys():
            with gen.Capturing() as output:
                gen.printing_class(expgraphchar[modexpn])
            print "'" + modexpn + "'; " + ', '.join(output)

    if debug: print "    End of '" + fname + "' "

    return modgraphchar, expgraphchar

def exp_headers(mod,config):
    """ Function to provide the headers and the experiments of a given model
      model= model
      config= configuration of the experiment
    """
    fname = 'exp_headers'

    # No CASE in python!
    #  case ${mod} in ...
    if mod == 'WRF':
        expers = config['WRFexps'].split(':')
        fhs = config['WRFheaders'].split(':')
    elif mod == 'LMDZ':
        expers = config['LMDZexps'].split(':')
        fhs = config['LMDZheaders'].split(':')
    elif mod == 'WRF_LMDZ':
        expers = config['WRF_LMDZexps'].split(':')
        fhs = config['WRF_LMDZheaders'].split(':')
    elif mod == 'DYNAMICO':
        expers = config['DYNAMICOexps'].split(':')
        fhs = config['DYNAMICOheaders'].split(':')
    elif mod == 'mthDYNAMICO':
        expers = config['mthDYNAMICOexps'].split(':')
        fhs = config['mthDYNAMICOheaders'].split(':')
    else:
        print errmsg
        print "  " + fname + ": model '" + mod + "' not ready!!"
        quit(-1)

    return expers, fhs

class VariableInf(object):
    """ Class which holds all information of a given variable
      name= CF name of the variabgle
      header= header of the file from which it can be computed
      model= file-header variable from which it can be directly computed
      diag= file-header diagnostics from which it can be computed
    """
    def __init__( self, name, fheader, model, diag):
       self.name= None
       self.fheader= None
       self.model= None
       self.diag= None
       if name is not None:
           self.name = name
           self.fheader= fheader
           self.model= model
           self.diag= diag

class ModelInf(object):
    """ Class which holds all information from a given model
      name= Acronym of the model
      model= Long description of the model
      dim[x/y/z/t/s]n= name of the dimensions of the model
      vdim[x/y/z/t/s]n= name of the variable-dimensions of the model
      dimensions= list of the dimensions
      vardimensions= list of the variable-dimensions
      timemodif= time modification to apply to the model files ('None' for nothing)
        'setorigin',[YYYYMMDDHHMISS]: re-set origin of times at [YYYYMMDDHHMISS]
    """
    def __init__( self, name, model, dx, dy, dz, dt, ds, vdx, vdy, vdz, vdt, vds,    \
      tmod):
        self.name= None
        self.model= None
        self.dimxn= None
        self.dimyn= None
        self.dimzn= None
        self.dimtn= None
        self.dimsn= None
        self.vardxn= None
        self.vardyn= None
        self.vardzn= None
        self.vardtn= None
        self.vardsn= None
        self.timemodif = None
        if name is not None:
            self.name = name
            self.model= model
            self.dimxn= dx
            self.dimyn= dy
            self.dimzn= dz
            self.dimtn= dt
            self.dimsn= ds
            self.vardxn= vdx
            self.vardyn= vdy
            self.vardzn= vdz
            self.vardtn= vdt
            self.vardsn= vds
            if tmod is not None:
                self.timemodif = tmod
            self.dimensions = [dx, dy, dz, dt, ds]
            self.vardimensions = [vdx, vdy, vdz, vdt, vds]

def variable_compute(idir,var,ftests,db):
    """ Function to retrieve the computation way of a given variable using a series of test files
      iwdir= directory with the test files
      var= name of the variable to test
      filtests= dictionary with the test files for each header
    """
    fname='variable_compute'

    cancompute = True
    for headerf in ftests.keys():
        filen=ftests[headerf]
        try:
            with gen.Capturing() as output:
                vmod, vdiag = ncvar.computevar_model(var, idir + '/' + filen)
        except:
            print errmsg
            print 'ncvar.computevar_model(' + var + ', ' + idir + '/' + filen + ')'
            for s1out in output: print s1out
            quit(-1)

        if vmod is None and vdiag is None:
            cancompute = False
        else:
            cancompute = True
            # Should be considered that variable can also be computed by both ways?
            break

    if not cancompute:
        print warnmsg
        print '  ' + fname + ": there is no way to compute '" + var +                \
          "' for model '" + mod
# Too extrict!
#        quit(-1)

    # Getting only the first possibility of model and diagnostic
    if vmod is not None:
        modV = vmod[0]
    else:
        modV = None
    if vdiag is not None:
        diagV = vdiag[0]
    else:
        diagV = None

    varcomp = VariableInf(var, headerf, modV, diagV)

    # a ';' list 'varcompute' it is created for each variable giving:
    #   [var]|[vark]|[headerf][varmod]|[vardiag]
    # This list will be used to compute a new file for each variable

    if db:
        print 'Variable information _______'
        gen.printing_class(varcomp)

    if db: print "    End of '" + fname + "' "

    return varcomp

def get_operations_var(vc,db):
    """ Function to provide the operations to make for each variable from 'VAR_' dictionary
      vc= 'VAR_' dictionary, dictionary with all the parameters started with 'VAR_'
    """
    fname = 'get_operations_var'

    LH = len('VAR_')
    iop = 1
    # list of operations
    doopers = {}
    # operations by variable
    operationsvar = {}
    # individual operations by variable (not repeated)
    indivoperationsvar = {}

    # Variables with a calculation whic requires vertical interpolation to all the 
    #   file (operations stating by 'VAR_pinterp')
    LVp = len('VAR_pinterp')
    varglobalp = []

    for oper in vc.keys():
        Loper = len(oper)
        opn = oper[LH:Loper+1]
        vns = vc[oper].split(':')
        ops1 = opn.split('+')
        doopers[opn] = ops1
        for vn in vns:
            if not operationsvar.has_key(vn):
                operationsvar[vn] = [opn]
                indivoperationsvar[vn] = ops1
            else:
                opers = operationsvar[vn]
                opers.append(opn)
                operationsvar[vn] = opers
                opers1 = indivoperationsvar[vn]
                for op1 in ops1:
                    if not gen.searchInlist(opers1, op1):
                        opers1.append(op1)
                indivoperationsvar[vn] = opers1

        if oper[0:LVp] == 'VAR_pinterp':
            for vn in vns:
                if not gen.searchInlist(varglobalp,vn): varglobalp.append(vn)

        iop = iop + 1

    if db:
        print '  operations to make _______'
        gen.printing_dictionary(doopers)
        print '  operations by variable _______'
        gen.printing_dictionary(operationsvar)
        print '  individual operations by variable _______'
        gen.printing_dictionary(indivoperationsvar)
        print '  variables with the vertical interpolation for all the file _______'
        print '    #', varglobalp

    if db: print "    End of '" + fname + "' "

    return doopers, operationsvar, indivoperationsvar, varglobalp

def pinterp_var(oper):
    """ Function to retrieve characteristics of the vertical interpolation for the operation
      oper= operation
        Wether vertical interpolation is:
          # 'global': for all file
          # 'local': at a given step of the process
          # 'none': no vertical interpolation
    >>> pinterp_var('last+pinterp+xmean')
    local
    >>> pinterp_var('pinterp+tmean+xmean')
    global
    >>> pinterp_var('xmean')
    none
    """
    fname = 'pinterp_var'

    pinn = 'pinterp'

    Lpinterp = len(pinn)
    if oper[0:Lpinterp] == pinn: 
        pkind = 'global'
        return pkind

    if oper.find(pinn) != -1:
        pkind = 'local'
    else:
        pkind = 'none'
        
    return pkind

def compvars_listconstruct(config, minf, Files, TestFiles, idir, odir, debug):
    """ Function to construct the list of variables to compute
      config= dictionary with the configuration of the execution
         variables to compute are get with all values started by `VAR_' in 'module_graphics.dat', and then
           providing a consecutive number of calculations separated by '+'
             VAR_[calc1]+[calc2] = tas:wss
           will compute first [calc1] and then [calc2] for 'tas' and 'wss'
      minf= class with information about the model
      Files= dictionary of files for heach header
      TestFiles= dictionary of files to test how to compute the variable
      odir= output directory
    """
    fname='compvars_listconstruct'

    varcomp = gen.get_specdictionary_HMT(config,H='VAR_')

    if debug:
        print '  variables to compute ________'
        gen.printing_dictionary(varcomp)

    # Getting operations by variable
    opers, varoper, indivaroper, varglobp = get_operations_var(varcomp,debug)    

    strfiles = gen.dictKeysVals_stringList(Files)
    strtestfiles = gen.dictKeysVals_stringList(TestFiles)

    Svarcompute = ''

    # Main dictionary with all the variables and their calculation
    allvarcomp = {}

    ivop = 0
    for vn in varoper:
        vcomp = variable_compute(idir,vn,TestFiles,False)
        for op in varoper[vn]:
            # Creation of a String as: [CFname]|[operation]|[fheader]|[vmodel]|[vdiag]|[globalP]
            #   [CFname]: CF-name of the variable (must appear in 'variables_values.dat')
            #   [operation]: operation to compute
            #   [fheader]: header of file with the required variables
            #   [vmodel]: Variable from model output
            #   [vdiag]: Varible as diagnostic from different variables of the model output
            #   [globalP]: Wether vertical interpolation is:
            #     'global': for all file
            #     'local': at a given step of the process
            #     'none': no vertical interpolation
            globalP = pinterp_var(op)
            if vcomp.model is not None:
                if type(vcomp.model) == type(list([1, 2])):
                    Smodel = ':'.join(vcomp.model)
                elif type(vcomp.model) == type('A'):
                    Smodel = vcomp.model
            else:
                Smodel = 'None'
            if vcomp.diag is not None:
                if type(vcomp.diag) == type(list([1, 2])):
                    Sdiag = ':'.join(vcomp.diag)
                elif type(vcomp.diag) == type('A'):
                    Sdiag = vcomp.diag
            else:
                Sdiag = 'None'

            # To get a combination of operations, all the precedent steps are kept
            #   thus, they are also available !
            seqops = op.split('+')
            for seqop in seqops:
                if seqop == seqops[0]: seqopS=seqop
                else: seqopS = seqopS + '+' + seqop 
                allvarcomp[vn+'_'+seqopS] = [vcomp.fheader, vcomp.model, vcomp.diag, \
                  globalP]

                Svc = vcomp.name + '|' + seqopS + '|' + vcomp.fheader + '|' +        \
                  Smodel + '|' + Sdiag + '|' + globalP
                if ivop == 0:
                    Svarcompute = Svc 
                else:
                    Svarcompute = Svarcompute + ',' + Svc 

                ivop = ivop + 1

    if debug:
        print '  Variables to compute _______'
        gen.printing_dictionary(allvarcomp)

    # Outwritting the varcompute to avoid next time (if it is not filescratch!)
    objf = open(odir + '/varcompute.inf', 'w')
    objf.write('files: ' + strfiles + '\n')
    objf.write('testfiles: ' + strtestfiles + '\n')
    objf.write('varcompute: ' + Svarcompute + '\n')
    objf.write('itotv: ' + str(ivop) + '\n')
    objf.close() 

    if debug: print "    End of '" + fname + "' "

    return allvarcomp, ivop

def read_varcomp_file(filen):
    """ Function to read 'varcompute.inf' and reconstruct the dictionaries
    """
    fname = 'read_varcomp_file'

    allvarcomp = {}
    
    objf = open(filen, 'r')

    for line in objf:
        vals = line.replace('\n','').split(' ')
        if vals[0] == 'files:':
            Files = gen.stringList_dictKeysVals(vals[1])
        elif vals[0] == 'testfiles:':
            TestFiles = gen.stringList_dictKeysVals(vals[1])
        elif vals[0] == 'varcompute:':
            for VvalsS in vals[1].split(','):
                
                Vvals = VvalsS.split('|')
                if Vvals[3] == 'None':
                    mod = None
                else:
                    mod = Vvals[3].split(':')
                if Vvals[4] == 'None':
                    diag = None
                else:
                    diag = Vvals[4].split(':')

                allvarcomp[Vvals[0]+'_'+Vvals[1]] = [Vvals[2], mod, diag, Vvals[5]]

#                # To get a combination of operations, all the precedent steps are kept
#                #   thus, they are also available !
#                seqops = Vvals[1].split('+')
#                for seqop in seqops:
#                    if seqop == seqops[0]: seqopS=seqop
#                    else: seqopS = seqopS + '+' + seqop 
#                    allvarcomp[Vvals[0]+'_'+seqopS] = [Vvals[2], mod, diag, Vvals[5]]
        elif vals[0] == 'itotv:':
            ivop = int(vals[1])

    objf.close()

    return Files, TestFiles, allvarcomp, ivop

def pinterpS(gP):
    """ For that variables which require vertical interpolation 'p' suffix to the file header is added
      gP= kind of pinterp: 'global', 'local', 'none'
    """
    fname = 'pinterpS'

    if gP != 'none':
        gPS = 'p'
    else:
        gPS = ''

    return gPS

def compute_variable(minf, idir, usefiles, odir, cvar, gP, scr, pyH, Tref, Tunits, db):
    """ Function to compute a variable
      minf= class with the information of the model
      idir= directory with the input files
      usefiles= dictionary of files as dict([headerf]) = [file1], ..., [fileN]
      odir= directory to write the output files
      cvar= class with the information of the variable: 'name', 'fheader', 'varmod', 'vardiag'
      gP= kind of vertical interpolation ('global', 'local', 'none')
      scr= should it be done from the scratch?
      pyH= location of the python HOME
      Tref= CF time reference
      Tunits= CF time units
    """
    fname='compute_variable'

    CFvarn=cvar.name
    headerf = cvar.fheader 
    modvar = cvar.model
    diagvar = cvar.diag

    cfiles = usefiles[headerf]

    # dimensions
    dnx = minf.dimxn
    dny = minf.dimyn
    # var-dimensions
    vdnx = minf.vardxn
    vdny = minf.vardyn

    # Computing separately and then joinging for all files
    Ntotfiles = len(cfiles)
    Nzeros = len(str(Ntotfiles))
    NStot = str(Ntotfiles).zfill(Nzeros)
    # For that variables which require vertical interpolation 'p' suffix to the 
    #   file header is added
    SgP = pinterpS(gP)

    # Getting in working dir
    os.chdir(odir)

    # File to keep track of all operations
    otrackf = open( odir + '/all_computevars.inf', 'a')

    # Computing variable
    ifile=1
    for cf in cfiles:
        ifS = str(ifile).zfill(Nzeros)
        ifilen = odir + '/' + CFvarn + '_' + headerf + SgP + '_' + ifS + '-' +       \
          NStot + '.nc'
        fileon = odir + '/' + CFvarn + '_' + headerf + SgP + '.nc'

        if scr:
            sout = sub.call('rm ' + ifilen + ' >& /dev/null', shell=True)
            sout = sub.call('rm ' + fileon + ' >& /dev/null', shell=True)

        if not os.path.isfile(ifilen) and not os.path.isfile(fileon):
# Since model direct values are retrieved from `variables_valules.dat' which was 
#   initially coincived as a way to only give variable attributes, range and color 
#   bars, if a variable has a diagnostic way to be computed, the later one will be 
#   preferred
            if db:
                print '  ' + fname + ": creation of variable file '" + CFvarn +      \
                  "' in file '" + ifilen + "' ..."
                print '    model variable:', modvar
                print '    diagostic variable:', diagvar

            if diagvar is None:
                # model variable
                if db: print '  '+fname+ ": variable direct from model '" +modvar+ "'"
                values = modvar + ',0,-1,-1'
                vs = modvar + ',' + vdnx + ',' + vdny + ',' + vdnz + ',' + vdnt

                try:
                    with gen.Capturing() as output:
                        ncvar.DataSetSection_multivars(values, idir + '/' + cf, vs)
                except:
                    print errmsg
                    print 'ncvar.DataSetSection_multivars('+values+', '+idir+'/'+cf+ \
                      ', '+vs+')'
                    for s1out in output: print s1out
                    quit(-1)

                newfile, loc = gen.search_sec_list(output,'succesfull')
                ofile = newfile[0].split(' ')[7]
                sub.call('mv ' + ofile + ' ' + ifilen, shell=True)

                # Keeping track of the operations
                pyins = pyH + "/nc_var.py -f " + idir + '/' + cf +                   \
                  " -o DataSetSection_multivars -v " + vs + " -S '" + values + "'"
                otrackf.write('\n')
                otrackf.write('# ' + CFvarn + " " + modvar + '\n')
                otrackf.write('python ' + pyins + '\n')

                # CF renaming of variable
                ncvar.chvarname(CFvarn,ifilen,modvar)
            else:
                # diagnostic variable
                if db: print '  '+fname+ ": variable as diagnostic '" +diagvar[0]+ "'"
                if minf.name == 'WRF':
                    dims = dnt+'@WRFtime,'+dnz+'@'+vdnz+','+dny+'@'+vdny+','+dnx+    \
                      '@'+vdnx
                else:
                    dims = dnt+'@'+vdnt+','+dnz+'@'+vdnz+','+dny+'@'+vdny+','+dnx+   \
                      '@'+vdnx
                diagn = diagvar[0]
                Ndiagvars = len(diagvar)
                diagc = '@'.join(diagvar[1:])
            
                values = '-f ' + idir + '/' + cf + " -d '" + dims + "' -v '" +       \
                  diagn + '|' + diagc + "'"
                try:
                    with gen.Capturing() as output:
                        sout = sub.call('python ' +pyH+ '/diagnostics.py ' + values, \
                          shell=True)
                except:
                    print errmsg
                    print 'python ' + pyH + '/diagnostics.py ' + values
                    for s1out in output: print s1out
                    quit(-1)
                if db:
                    for s1out in output: print s1out

                sout = sub.call(' mv diagnostics.nc ' + ifilen, shell=True)

                # Keeping track of the operations
                pyins = 'python ' + pyH + '/diagnostics.py ' + values
                otrackf.write('\n')
                otrackf.write('# ' + CFvarn + " " + diagn + '\n')
                otrackf.write(pyins + '\n')

            # Attaching necessary variables for the pressure interpolation
            if gP != 'none':
                if minf.name == 'WRF' or minf.name == 'WRF_LMDZ':
                    requiredinterpvars = ['P', 'PB', 'PSFC', 'PH', 'PHB', 'HGT', 'T',\
                      'QVAPOR', 'XLONG', 'XLAT', 'Times']
                elif minf.name == 'LMDZ':
                    requiredinterpvars = ['pres', 'psol', 'geop', 'phis', 'temp',    \
                      'ovap', 'lon', 'lat', 'time_counter']
                elif minf.name == 'DYNAMICO':
                    requiredinterpvars = ['pres', 'psol', 'geop', 'phis', 'temp',    \
                      'ovap', 'lon', 'lat', 'time_counter']
                elif minf.name == 'mthDYNAMICO':
                    requiredinterpvars = ['pres', 'psol', 'geop', 'phis', 'temp',    \
                      'ovap', 'lon', 'lat', 'time_centered']
                else:
                    print errmsg
                    print fname + ": for model '" + minf.name + "' required " +      \
                      " variables for vertical interpolation are not known !!"
                    quit(-1)

                print "   " + fname + ": adding variables:", requiredinterpvars,     \
                  ' to allow pressure interpolation'
                for rqv in requiredinterpvars:
                    try:
                        with gen.Capturing() as output:
                            ncvar.fvaradd(idir+'/'+cf+','+rqv,ifilen)
                    except:
                        print errmsg
                        print 'fvaradd('+idir+'/'+cf+','+rqv+', '+ifilen+')'
                        for s1out in output: print s1out
                        quit(-1)
                    if db:
                        for s1out in output: print s1out

            # CFification of files
            if minf.name == 'WRF' or minf.name == 'WRF_LMDZ':
                # Assuming that all diagnostics are outputed with CF-times
                values = vdnx + ':' + vdny + ':' + 'time:'+ Tref + ':' + Tunits
                try:
                    with gen.Capturing() as output:
                        ncvar.WRF_to_newCF(values, ifilen, 'all') 
                        sub.call('mv CF_WRF.nc '+ifilen + ' >& /dev/null', shell=True)
                except:
                    print errmsg
                    print 'WRF_to_newCF('+ values +', ' + ifilen + ',all)'
                    for s1out in output: print s1out
                    # Removing file in order to make sure that it will be redone
                    print '  ' + fname + " removing file '" + ifilen + "' ..."
                    sub.call('rm ' + ifilen + ' >& /dev/null', shell=True)
                    quit(-1)
                if db:
                    for s1out in output: print s1out
            elif minf.name == 'LMDZ':
                try:
                    with gen.Capturing() as output:
                        ncvar.LMDZ_toCF(Tunits + '!since!' + Tref, ifilen)
                except:
                    print errmsg
                    print 'LMDZ_toCF('+ Tunits +'!since!'+ Tref + ', ' + ifilen + ')'
                    for s1out in output: print s1out
                    # Removing file in order to make sure that it will be redone
                    print '  ' + fname + " removing file '" + ifilen + "' ..."
                    sub.call('rm ' + ifilen + ' >& /dev/null', shell=True)
                    quit(-1)
                if db:
                    for s1out in output: print s1out
            # At this stage, let's be simple. It might be necessary to complicate it! Since
            #   DYNAMICO might have 'time_counter' or 'time_instant'
            elif minf.name == 'DYNAMICO':
                try:
                    with gen.Capturing() as output:
                        ncvar.DYNAMICO_toCF(Tunits + '!since!' + Tref, ifilen)
                except:
                    print errmsg
                    print 'DYNAMICO_toCF('+Tunits+'!since!'+Tref+ ', ' + ifilen + ')'
                    for s1out in output: print s1out
                    # Removing file in order to make sure that it will be redone
                    print '  ' + fname + " removing file '" + ifilen + "' ..."
                    sub.call('rm ' + ifilen + ' >& /dev/null', shell=True)
                    quit(-1)
                if db:
                    for s1out in output: print s1out
            elif minf.name == 'mthDYNAMICO':
                try:
                    with gen.Capturing() as output:
                        ncvar.mthDYNAMICO_toCF(Tunits + '!since!' + Tref, ifilen)
                except:
                    print errmsg
                    print 'mthDYNAMICO_toCF('+Tunits+'!since!'+Tref+ ', ' + ifilen + ')'
                    for s1out in output: print s1out
                    # Removing file in order to make sure that it will be redone
                    print '  ' + fname + " removing file '" + ifilen + "' ..."
                    sub.call('rm ' + ifilen + ' >& /dev/null', shell=True)
                    quit(-1)
                if db:
                    for s1out in output: print s1out
            else:
                print errmsg
                print '  ' + fname + ": no CFification for model '" +minf.name+ "' !!"
                print "    available ones: 'WRF', 'WRF_LMDZ', 'LMDZ', 'DYNAMICO'" +  \
                  ", 'mthDYNAMICO'"
                quit(-1)

        ifile = ifile + 1

    otrackf.close()

    # Joining variable files
    if not os.path.isfile(fileon):
        if Ntotfiles > 1:
            if db:
                print '  ' + fname + ": concatenating all variable files to " +      \
                  "create '" + fileon + "..."
            try:
                with gen.Capturing() as output:
                    ncvar.netcdf_fold_concatenation_HMT(odir+',time,time', CFvarn +  \
                      '_' + headerf + SgP +'_,-,.nc', 'all')
            except:
                print errmsg
                print 'netcdf_fold_concatenation_HMT('+odir+',time,time ' + CFvarn + \
                   '_' + headerf +SgP + '_,-,.nc, all)'
                for s1out in output: print s1out
                quit(-1)
            if db:
                for s1out in output: print s1out

            sout = sub.call('mv netcdf_fold_concatenated_HMT.nc '+fileon, shell=True)
            if os.path.isfile(fileon):
                sout = sub.call('rm '+CFvarn + '_' + headerf + SgP + '_*-*.nc >& ' + \
                  '/dev/null', shell=True)
        else:
            sout = sub.call('mv ' + CFvarn + '_' + headerf + SgP + '_1-1.nc ' +      \
                  fileon, shell=True)

    # Re-setting time of final concatenated file
    if minf.timemodif is not None:
        if db:
            print '  ' + fname + ": Modifying times of the file by '" +              \
              minf.timemodif + "' !!"
        try:
            with gen.Capturing() as output:
                ncvar.time_reset(minf.timemodif, fileon,'time')
        except:
            print errmsg
            print 'time_reset(' + minf.timemodif + ', ' + fileon + ', time)'
            for s1out in output: print s1out
            # Removing file in order to make sure that it will be redone
            print '  ' + fname + " removing file '" + fileon + "' ..."
            sub.call('rm ' + fileon + ' >& /dev/null', shell=True)
            quit(-1)
        if db:
            for s1out in output: print s1out

    if db: print "    End of '" + fname + "' "

    return

def compute_statistics(minf, config, idir, usefiles, odir, cvar, gP, Opers, scr, db):
    """ Function to compute different statistics it will take previous steps if they 
        are availale
      minf= class with the information of the model
      config= dictionary with the configuration of the experiment
      idir= directory with the input files
      usefiles= ',' list of files to use [file1],...,[fileN]
      odir= directory to write the output files
      cvar= class with the information of the variable: 'name', 'fheader', 'varmod', 'vardiag'
      gP= kind of vertical interpolation ('global', 'local', 'none')
      Opers= kind of operation: (as possible multiple consecutive combination of operations separated by '+'
        [calc1]+[calc2] will compute first [calc1] and then [calc2] 
        acc: temporal accumulated values
        diff: differences between models
        direct: no statistics
        last: last temporal value
        Lmean: latitudinal mean values
        Lsec: latitudinal section (latitudinal value must be given, [var]@[lat]) 
        lmean: longitudinal mean values
        lsec: longitudinal section (longitudinal value must be given, [var]@[lat])
        pinterp: pressure interpolation (to the given $plevels)
        smean: spatial-weighted mean values
        tmean: temporal mean values
        tstd: temporal standard deviation values
        tturb: Taylor's turbulence decomposition value (x - <x>) for time
        tvar: temporal variance values
        xmean: x-axis mean values
        xvar: x-axis variance values
        ymean: y-axis mean values
        zsum: vertical aggregated values
      scr= should it be done from the scratch?
    """
    fname='compute_statistics'

    # Getting a previous file name to continue with(for combinations)
    CFvarn=cvar.name
    headerf = cvar.fheader 
    modvar = cvar.model
    diagvar = cvar.diag

    cfiles = usefiles[headerf]

    # dimensions
    dnx = minf.dimxn
    dny = minf.dimyn
    # var-dimensions
    vdnx = minf.vardxn
    vdny = minf.vardyn
    # Model dimension-variable dictionary
    Modeldimvardict = {dnx: vdnx, dny: vdny, dnz: vdnz, dnt: vdnt}

    # Some experiment configuration values
    #  plevels= ':' separated list of pressures (in Pa) to use for the vertical interpolation
    plevels = config['plevels']
    #  pyH= location of the python HOME
    pyH = config['pyHOME']
    #  Tref= CF time reference
    Tref = config['CFreftime']
    #  Tunits= CF time units
    Tunits = config['CFunitstime']
    #  opsur = operation surnames
    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
    Opsurs = {}
    for opk in opsur:
        opn = opk.split('_')[1]
        vals = opsur[opk]
        Opsurs[opn] = vals.split(':')
#    if db:
#        print '  ' + fname + ' operation surnames _______'
#        gen.printing_dictionary(Opsurs)

    # Getting in working dir
    os.chdir(odir)

    # File to keep track of all operations
    otrackf = open( odir + '/all_statsvars.inf', 'a')

    # For that variables which require vertical interpolation 'p' suffix to the file
    #   header is added. 
    # After some operations some CF dimension-variables might not be kept

    if gP != 'none':
        SgP = 'p'
        varnCFs = ['lon', 'lat', 'pres', 'time']
    else:
        SgP = ''
        varnCFs = ['lon', 'lat', 'time']

    # CF dimension-variable dictionary
    CFdimvardict = {'lon': 'lon', 'lat': 'lat', 'pres': 'pres', 'time': 'time'}

    # Input file
    ifilen = odir + '/' + CFvarn + '_' + headerf + SgP + '.nc'

    # Computing variable statisitcs
    istats=0

    # List of not computed statistics
    wrongstats = []

    opers = Opers.split('+')
    Fopers = ''
    if db: print "    computing statistics of variable:'", CFvarn, "' operation '" + \
      Opers + "' using file '" + ifilen + "...'"
    for op in opers:

        # File name from previous operations, and name of the variable within the 
        #   file (some operations change the name of the variable)
        if op == opers[0]:
            Fopers = op
            prevfile = ifilen
            vninF = CFvarn
        else:
            Fopers = Fopers + '_' + op
            prevfile = fileon
        fileon = odir + '/' + CFvarn + '_' + headerf + SgP + '_' + Fopers + '.nc'

        # Adding required variables for the vertical interpolation in all operations
        SvarnCFs = ',' + ','.join(varnCFs)

        if gP != 'none':
            if gP == 'local':
                CFvarnp = vninF + ',P,PB,PSFC,PH,PHB,HGT,T,QVAPOR,XLONG,XLAT,Times'+ \
                  SvarnCFs
            else:
                CFvarnp = vninF + SvarnCFs
        else:
            CFvarnp = vninF + SvarnCFs

        if scr:
            sout = sub.call('rm ' + fileon + ' >& /dev/null', shell=True)

        hprevfile = prevfile[0:len(prevfile)-3]

        if db: print '      ', op, "using '" + prevfile + "' "
        # Just produce file in case output file does not exist
        if not os.path.isfile(fileon):
            if db: 
                print '  ' + fname + ": creation of file '" + fileon + "' ..."
            dofile = True
        else:
            dofile = False
        
        # Variables to be kept in the final file
        varkeep = []

        if op == 'acc':
            # temporal accumulated values
            print "  " + fname + ": kind '" + op + "' not ready !!"
            wrongstats.append(CFvarn + '_' + opers)
            break

        elif op[0:11] == 'dimvarvalue':
            # Slicing along the index at the nearest given value of a dimension-variable
            dimvarn = op.split('~')[1]
            dimvalue = op.split('~')[2]
            vals = dimvarn + ',' + dimvalue + ',' + dimvalue + ',0'
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.DataSetSection_multivars(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'DataSetSection_multivars('+vals+',', prevfile, ',' +      \
                      CFvarnp + ')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv DataSetSection_multivars.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o DataSetSection_multivars -S '" + vals +   \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + CFvarnp + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing the given variable-dimension
            varnCFs.remove(dimvarn)

        elif op == 'direct':
            # no statistics
            if dofile:
                sout = sub.call('mv ' + prevfile + ' ' + fileon, shell=True)
        elif op == 'last':
            # last temporal value
            vals='time,-9,0,0'
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.DataSetSection(vals,prevfile)
                except:
                    print errmsg
                    print 'DataSetSection('+vals+',', prevfile, ')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out
                
                sout = sub.call('mv ' + hprevfile + '_time_B-9-E0-I0.nc ' + fileon, \
                  shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o DataSetSection -S '" + vals + "' -f " +  \
                  prevfile
                otrackf.write("\n")
                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

        elif op =='Lmean':
            # latitudinal mean values
            print "  " + fname + ": kind '" + op + "' not ready !!"
            wrongstats.append(CFvarn + '_' + opers)
            break
        elif op =='Lsec':
            # latitudinal section (latitudinal value must be given, [var]@[lat]) 
            print "  " + fname + ": kind '" + op + "' not ready !!"
            wrongstats.append(CFvarn + '_' + opers)
            break
        elif op =='lmean':
            # longitudinal mean values
            print "  " + fname + ": kind '" + op + "' not ready !!"
            wrongstats.append(CFvarn + '_' + opers)
            break
        elif op =='lsec':
            # longitudinal section (longitudinal value must be given, [var]@[lon])
            print "  " + fname + ": kind '" + op + "' not ready !!"
            wrongstats.append(CFvarn + '_' + opers)
            break
        elif op == 'pinterp':
            # pinterp: pressure interpolation (to the given $plevels)
            vals=plevels + ',1,1'
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.pinterp(vals,prevfile,vninF)
                except:
                    print errmsg
                    print 'pinterp('+vals+',', prevfile, ','+vninF+')'
                    for s1out in output: print s1out
                    ncvar.pinterp(vals,prevfile,vninF)
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv pinterp.nc ' + fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o pinterp -S '" + vals + "' -f " +          \
                  prevfile + "-v " + vninF
                otrackf.write("\n")
                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

                # adding CF lon,lat,time in WRF files
                if minf.name == 'WRF' or minf.name == 'WRF_LMDZ':
                    values = vdnx + ':' + vdny + ':'+ Tref + ':' + Tunits
                    if db:
                        print "   CFification of '" + fileon + "' with:", values, "..."
                    try:
                        with gen.Capturing() as output:
                            ncvar.WRF_toCF(values, fileon)
                    except:
                        print errmsg
                        print 'ncvar.WRF_toCF(' + values+ ', ' + fileon + ')'
                        for s1out in output: print s1out
                        # Removing file in order to make sure that it will be redone
                        print '  ' + fname + " removing file '" + fileon + "' ..."
                        sub.call('rm ' + fileon + ' >& /dev/null', shell=True)
                        quit(-1)

            # vertical interpolation variables are no more needed
            CFvarnp = vninF
            gP = 'none'

        elif op == 'smean':
            # spatial-weighted mean values
            noLlvdims = list(varnCFs)
            noLlvdims.remove('lon')
            noLlvdims.remove('lat')
            Saddvars = ':'.join(noLlvdims)
            if len(Saddvars) > 1:
                vals = 'reglonlat,lon,lat,lon,lat,' + Saddvars
            else:
                vals = 'reglonlat,lon,lat,lon,lat,None'

            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.SpatialWeightedMean(vals,prevfile,vninF)
                except:
                    print errmsg
                    print 'SpatialWeightedMean('+vals+',', prevfile, ','+vninF+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                # renaming variable
                ncvar.chvarname(vninF+'mean', 'SpatialWeightedMean_reglonlat.nc',    \
                  vninF+'spaceweightmean')
                sout = sub.call('mv SpatialWeightedMean_reglonlat.nc '+fileon,       \
                  shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o SpatialWeightedMean -S '" + vals +        \
                  "' -f " + prevfile + " -v " + vninF
                otrackf.write("\n")
                otrackf.write("# " + vninF + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'lon', 'lat'
            varnCFs.remove('lon')
            varnCFs.remove('lat')

            varkeep.append(vninF + 'mean')

        elif op == 'tmean':
            # temporal mean values
            vals='time|-1,time,mean,' + ':'.join(varnCFs) + ':' + vdnz
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_mean.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'time'
            varnCFs.remove('time')

            varkeep.append('timestats')

        elif op == 'tstd':
            # temporal standard deviation values
            vals='time|-1,time,std,' + ':'.join(varnCFs) + ':' + vdnz
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_std.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'time'
            varnCFs.remove('time')

            varkeep.append('timestats')

        elif op == 'tturb':
            # temporal turbulent values
            vals='time|-1,time,turb,' + ':'.join(varnCFs) + ':' + vdnz
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_turb.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'time'
            varkeep.append('timestats')

        elif op == 'tvar':
            # temporal variance values (Taylor's turbulence)
            vals='time|-1,time,var,' + ':'.join(varnCFs) + ':' + vdnz
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_var.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'time'
            varnCFs.remove('time')

            varkeep.append('timestats')

        elif op == 'xmean':
            # x-axis mean values
            vals='lon|-1,lon,mean,' + ':'.join(varnCFs) + ':' + vdnz
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_mean.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'lon'
            varnCFs.remove('lon')

            varkeep.append('lonstats')

        elif op == 'xvar':
            # x-axis variance values
            vals='lon|-1,lon,var,' + ':'.join(varnCFs) + ':' + vdnz
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_var.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'lon'
            varnCFs.remove('lon')

            varkeep.append('lonstats')

        elif op == 'ymean':
            # y-axis mean values
            vals='lat|-1,lat,mean,' + ':'.join(varnCFs) + ':' + vdnz
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_mean.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + CFvarn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'lat'
            varnCFs.remove('lat')

            varkeep.append('latstats')

        elif op == 'zsum':
            # vertical aggregated values
            print "  " + fname + ": kind '" + op + "' not ready !!"
            wrongstats.append(CFvarn + '_' + opers)
            break
        else:
            print errmsg
            print '  ' + fname + ": operation '" + op + "' not ready !!"
            quit(-1)

        # End of kind of operation

        # Variable name in file (vninF) changed due to operation 
        #   but only if previous operation does not the same 'statistic'
        chvn = gen.dictionary_key_list(Opsurs, op)
        if chvn is not None:
            oldvninF = vninF
            vninF = vninF + chvn
            CFvarnp = CFvarnp.replace(oldvninF,vninF)
                
        if dofile:
            if len(varkeep) > 0:
                varkeepS = ',' + ','.join(varkeep)
            else:
                varkeepS = ''

            # Adding such CF dimension variables which might not be in the 
            #   post-operation file
            try:
                with gen.Capturing() as output:  
                    varinfile = ncvar.ivars(fileon)
            except:
                print errmsg
                print 'ivars('+fileon+')'
                for s1out in output: print s1out
                # Removing file in order to make sure that it will be redone
                print '  ' + fname + " removing file '" + fileon + "' ..."
                sub.call('rm ' + fileon + ' >& /dev/null', shell=True)
                quit(-1)

            for CFvn in varnCFs:
                if not gen.searchInlist(varinfile, CFvn):
                    if db:
                        print '  ' + fname + ": recupering CF variable '" + CFvn + "'"
                    try:
                        with gen.Capturing() as output:
                            ncvar.fvaradd(prevfile+','+CFvn, fileon)
                    except:
                        print errmsg
                        print 'fvaradd(' + prevfile + ',' + CFvn +', '+fileon+')'
                        ncvar.fvaradd(prevfile+','+CFvn, fileon)
                        for s1out in output: print s1out
                        # Removing file in order to make sure that it will be redone
                        print '  ' + fname + " removing file '" + fileon + "' ..."
                        #sub.call('rm ' + fileon + ' >& /dev/null', shell=True)
                        quit(-1)

                    if db:
                        for s1out in output: print s1out

            totalvarkeeps = CFvarnp + ',' + ','.join(varnCFs) + varkeepS
            try:
                with gen.Capturing() as output:
                    oclean = ncvar.cleaning_varsfile(totalvarkeeps,fileon)
            except:
                print errmsg
                print 'cleaning_varsfile('+totalvarkeeps+','+fileon+')'
                ncvar.cleaning_varsfile(totalvarkeeps,fileon)
                for s1out in output: print s1out
                # Removing file in order to make sure that it will be redone
                print '  ' + fname + " removing file '" + fileon + "' ..."
                sub.call('rm ' + fileon + ' >& /dev/null', shell=True)
                quit(-1)

            if db:
                for s1out in output: print s1out

    # End of operations
    if db:
        print fname + ": successful calculation of '" + vninF + "' from '" +        \
          CFvarn + "' '" + Opers + "' !!"

    if len(wrongstats) > 1:
        print warnmsg
        print '  ' + fname + ": statisitcs not possible to compute:", wrongstats

    if db: print "    End of '" + fname + "' "

    return 

def compute_vars(config, modinf, idir, odir, Files, allvarcomp, fscr, debug):
    """ Function to compute the variables
      config= Configuration of the experiment
      modinf= class with information about the model
      idir= input experiment folder
      odir= output experiment folder
      Files= dictionary with the header and the files' correspondence
      allvarcomp= dictionary with all the variables to compute and their information
      fscr= whether files should be done from the scratch or not
    """
    fname = 'compute_vars'

    for vopn in allvarcomp.keys():
        # variable & operation
        vn = vopn.split('_')[0]
        op = vopn.split('_')[1]
        vopnV = allvarcomp[vopn]
        fheader = vopnV[0]
        if type(vopnV[1]) == type([1, 2]) and len(vopnV[1]) == 1:
            vals = vopnV[1]
            model = vals[0]
        else:
            model = vopnV[1]
        diag = vopnV[2]
        globalP = vopnV[3]

        if debug:
            print '      ' + vn + ': ' + op

        # Computing CF variable
        if model is not None or diag is not None:
            vinf = VariableInf(vn,fheader,model,diag)
            if op.find('pinterp') != -1: 
                SP = 'p'
            else:
                SP = ''

            # Comppute variable
            finalfilen = odir + '/' + vn + '_' + fheader + SP + '.nc'
            if fscr:
                sout = sub.call('rm ' + finalfilen + ' >& /dev/null', shell=True)

            if not os.path.isfile(finalfilen):
                compute_variable(modinf, idir, Files, odir, vinf, globalP, fscr,     \
                  config['pyHOME'], config['CFreftime'], config['CFunitstime'], debug)

            # Compute variable statistics
            finalfilen = odir + '/' + vn + '_' + fheader + SP + '_' +                \
              op.replace('+','_') + '.nc'
            if fscr:
                sout = sub.call('rm ' + finalfilen + ' >& /dev/null', shell=True)

            if not os.path.isfile(finalfilen):
                compute_statistics(modinf, cnf, iwdir, Files, owdir, vinf, globalP,  \
                  op, fscr, debug)

        else:
            print errmsg
            print '  ' + fname + ": neither 'model' or 'diag' variables for '" + vn  \
              + "' !!"
            quit(-1)

    return

def get_plots_var(vc,pltk,db):
    """ Function to provide the plots to make for each variable from pltk ('DIRPLT_', or 'DIFFPLT_') dictionary
      vc= pltk dictionary, dictionary with all the parameters started with pltk
      pltk= kind of plots. Values started by 'DIRPLT_', or 'DIFFPLT_' in model configuration
    """
    fname = 'get_plots_var'

    LH = len(pltk)
    # list of plots
    doplots = {}
    # plots by variable
    plotsvar = {}
    # individual plots by variable (not repeated)
    indivplotsvar = {}

    # Variables with a plot which requires vertical interpolation
    varplotp = []

    ipl = 0
    for plot in vc.keys():
        Lplot = len(plot)
        pln = plot[LH:Lplot+1]
        varops = vc[plot].split(':')
        for vn in varops:
            VOn = vn.split('#')
            # Variables and operations in plot
            vplt = []
            oplt = []
            voplt = []
            for vop in VOn:
                vv = vop.split('|')[0]
                oo = vop.split('|')[1]
                vplt.append(vv)
                oplt.append(oo)
                if not gen.searchInlist(voplt,vop): voplt.append(vop)
            vpltS = '#'.join(vplt)
            opltS = '#'.join(oplt)
            
            if not doplots.has_key(pln):
                doplots[pln] = [vn]
            else:
                plots = doplots[pln]
                plots.append(vn)
                doplots[pln] = plots

            if not plotsvar.has_key(vn):
                plotsvar[vn] = [pln]
            else:
                plots = plotsvar[vn]
                plots.append(pln)
                plotsvar[vn] = plots

            for VOv in voplt:
                if not indivplotsvar.has_key(VOv):
                    indivplotsvar[VOv] = [pln]
                else:
                    plots = indivplotsvar[VOv]
                    plots.append(pln)
                    indivplotsvar[VOv] = plots

        if vc[plot].find('pinterp') != -1:
            if not gen.searchInlist(varplotp,vn): varplotp.append(pln)

        ipl = ipl + 1

    if db:
        print '  plots to make _______'
        gen.printing_dictionary(doplots)
        print '  plots by variable|operation _______'
        gen.printing_dictionary(plotsvar)
        print '  individual plots by variable|operation _______'
        gen.printing_dictionary(indivplotsvar)
        print '  plots with the vertical interpolation _______'
        print '    #', varplotp

    if db: print "    End of '" + fname + "' "

    return doplots, plotsvar, indivplotsvar, varplotp

def plots_listconstruct(config, pkind, finf, odir, debug):
    """ Function to create the list of plots to draw
      config= Configuration of the experiment
      pkind= kind of plots ('DIRPLT', 'DIFFPLT')
      finf= file with the instructions of the plots 
      odir= output experiment folder
    """
    fname='plots_listconstruct'
  
    dirplot = gen.get_specdictionary_HMT(config, H=pkind+'_')

    # No plots because there are not plots of this kind
    if dirplot is None: 
        print warnmsg
        print '  ' + fname + ": no plots for type '" + pkind + "' !!"
        return {}, 0

    if debug:
        print "  '" + pkind + "' plots to draw ________"
        gen.printing_dictionary(dirplot)

    # Getting plots by variable
    plots, varplot, indivarplot, plotp = get_plots_var(dirplot, pkind+'_', debug)

    Svarplot = gen.dictKeysVals_stringList(plots,cV=':')

    Nplots = 0
    for pl in Svarplot.split(','):
        Nplots = Nplots + len(pl.split(':'))

    # Outwritting the varcompute to avoid next time (if it is not filescratch!)
    objf = open(finf, 'w')
    objf.write('plots: ' + Svarplot + '\n')
    objf.write('itotp: ' + str(Nplots) + '\n')
    objf.close() 

    if debug: print "    End of '" + fname + "' "

    return plots, Nplots

def read_plot_file(figfile):
    """ Function to read the file with the information about the  plots to draw
      figfile= file with the information
    """
    fname = 'read_plot_file'

    if not os.path.isfile(figfile):
        print errormsg
        print '  ' + fname + ": figures file '" + figfile + "' does not exist !!"
        quit(-1)

    objf = open(figfile, 'r')
    for line in objf:
        if line[0:1] != '#' and len(line) > 1:
            values = line.replace('\\n','').split(' ')
            if values[0] == 'plots:':
                if values[1] != 'none':
                    plots = gen.stringList_dictKeysVals(values[1],cV=':')
                else:
                    pots = {}
            elif values[0] == 'itotp:':
                Nplots = int(values[1])

    objf.close()

    return plots, Nplots

def varnoper(vn, oper, opsurdict):
    """ Function to provide name of the variable after operation as result of addition of the 'surnames'
      vn= variable name
      oper= '+' separated list of operations
      opsurdict= dictionary with the surname of operations
    >>> OpSurDict = {'mean': ['tmean', 'xmean', 'ymean']}
    >>> varnoper('tas', 'pinterp+tmean+xmean', OpSurDict)
    tasmeanmean
    """
    fname = 'varnoper'

    newvn = vn

    for op in oper.split('+'):
        surname = gen.dictionary_key_list(opsurdict,op)
        if surname is not None:
            newvn = newvn + surname

    return newvn

def draw_plot(kplot, CFvplot, fplot, vplot, dplot, pplot, finame, tfig, kfig, mapval,\
  tvals, expgraphc, od, pyH, fscr, db):
    """ Function to draw a plot
      kplot= kind of plot
      CFvplot= list of the CFnames of the variables
      fplot= list with the files in the plot
      vplot= list with the name of the variables in each file
      dplot= list of the dimensions in each file
      pplot= list of pictoric characteristics for each variable
      finame= name of the figure
      tfig= title of the figure
      kfig= kind of the figure
      mapval= value of the map
      tvals= list with time values
        tunits: units of time ('!' for spaces)
        timekind: kind of time ticks in the plot
        timefmt: format of time ticks in the plot
        timelabel: label of time-axis in the plot
      expgraphc= dictionary with expriment graphical characteristics
      od= output directory
      pyH= python HOME
      fscr= whether figure should be done from scratch
    """
    fname = 'draw_plot'

    if fscr:
        sout = sub.call('rm ' + finame + ' >& /dev/null', shell=True)

    os.chdir(od)

    # CF variable-dimensions
    CFvardims = {'lon':'lon', 'lat': 'lat', 'pres': 'pres', 'time': 'time'}

    # Tracking file with all the figures
    trkobjf = open('all_figures.inf', 'a')

    # Time-characteristics
    tunits = tvals[0]
    timekind = tvals[1]
    timefmt = tvals[2]
    timelabel = tvals[3]

    # Graphical labels and configuration
    

    if not os.path.isfile(finame):
        if db:
            print "    drawing '" + finame + "' ..."
            print '      plot kind:', kplot
            print '      CFvars:', CFvplot
            print '      files:', fplot
            print '      in file vars:', vplot
            print '      dims:', dplot
            print '      pictoric values:', pplot
            print '      figure name:', finame
            print '      title:', tfig.replace('!',' ')
            print '      kind:', kfig
            print '      map:', mapval
            print '      time units:', tunits

        if kplot == 'diffmap2Dsfc':
            print "  " + fname + ": kind of plot '" + kplot + "' not ready !!"
            quit(-1)
        elif kplot == 'diffmap2Dz':
            print "  " + fname + ": kind of plot '" + kplot + "' not ready !!"
            quit(-1)
        elif kplot == 'map2Dsfc':
            print "  " + fname + ": kind of plot '" + kplot + "' not ready !!"
            quit(-1)
        elif kplot == 'map3D':
            print "  " + fname + ": kind of plot '" + kplot + "' not ready !!"
            quit(-1)
        elif kplot == '2lines':
            figtit = tfig.replace('!','|')

            lineAstdn = CFvplot[0]
            lineBstdn = CFvplot[1]
            figfs = ','.join(fplot)

            lineA = pplot[0]
            Arange = str(lineA[0]) + ',' + str(lineA[1])
            # Values are changed from specific values in `model_graphics.dat'
            if lineA[3].find('%') == -1:
                Aline = lineA[2]
                Akind = lineA[3]
                Amark = lineA[4]
                Asize = lineA[5]
            else:
                Aline = 'red'
                Akind = '-'
                Amark = ','
                Asize = 2.

            lineB = pplot[1]
            Brange = str(lineB[0]) + ',' + str(lineB[1])
            # Values are changed from specific values in `model_graphics.dat'
            if lineB[3].find('%') == -1:
                Bline = lineB[2]
                Bkind = lineB[3]
                Bmark = lineB[4]
                Bsize = lineB[5]
            else:
                Bline = 'blue'
                Bkind = '-'
                Bmark = ','
                Bsize = '2.'

            # It is assumed that if the space variable is 'lon': is desired a 
            #   (lon, vals) plot if it is 'lat': then (vals, lat) plot
            if gen.searchInlist(dplot[0],'lon'):
                spacedim ='lon'
                axisvals = 'y'
                figmid = 'longitudinal'
            elif gen.searchInlist(dplot[0],'lat'):
                spacedim = 'lat'
                axisvals = 'x'
                figmid = 'meridional'
            elif gen.searchInlist(dplot[0],'pres'):
                spacedim = 'pres'
                axisvals = 'x'
                figmid = 'vertical'
            else:
                print errmsg
                print '  ' + fname + ": in '2lines' only ready for: 'lon', 'lat'," + \
                 " 'pres' as common dimension !!"
                print '  dimensions in plot:', dplot[0]
                quit(-1)

            figtit = '!'.join(tfig.split('!')[0:2]) + '!' + figmid + '!' +           \
              '!'.join(tfig.split('!')[2:])

            graphvals = spacedim + ':' + Arange + ':' + Brange + ':Extrs:' +         \
              axisvals + ':' + ','.join(CFvplot) + ':' + Aline + ',' + Bline +       \
              ':2.,2.:' + Akind + ',' + Bkind + ':2.,2.:,:' + figtit + ':' +         \
              spacedim + ':0|auto:' +  finame.replace('.'+kfig, '') + ':' + kfig +   \
              ':True'

            fvarS = ','.join(vplot)

            drwins = 'draw_2lines'
            plotins = 'python ' + pyH + '/drawing.py -f ' + figfs +' -o ' + drwins + \
              " -S '" + graphvals + "' -v " + fvarS

            try:
                with gen.Capturing() as output:
                    sout = sub.call(plotins, shell=True)
            except:
                print errmsg
                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
                print sout
                for s1out in output: print s1out
                quit(-1)

            # keeping all figures
            trkobjf.write('\n')
            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
            trkobjf.write(plotins + '\n')

        elif kplot == '2linesTime':
            figtit = tfig.replace('!','|')

            lineAstdn = CFvplot[0]
            lineBstdn = CFvplot[1]
            figfs = ','.join(fplot)

            lineA = pplot[0]
            Arange = str(lineA[0]) + ',' + str(lineA[1])
            # Values are changed from specific values in `model_graphics.dat'
            if lineA[3].find('%') == -1:
                Aline = lineA[2]
                Akind = lineA[3]
                Amark = lineA[4]
                Asize = lineA[5]
            else:
                Aline = 'red'
                Akind = '-'
                Amark = ','
                Asize = '2.'

            lineB = pplot[1]
            Brange = str(lineB[0]) + ',' + str(lineB[1])
            # Values are changed from specific values in `model_graphics.dat'
            if lineB[3].find('%') == -1:
                Bline = lineB[2]
                Bkind = lineB[3]
                Bmark = lineB[4]
                Bsize = lineB[5]
            else:
                Bline = 'blue'
                Bkind = '-'
                Bmark = ','
                Bsize = '2.'

            timeaxis = 'x'
            figmid = 'evolution'

            figtit = '!'.join(tfig.split('!')[0:2]) + '!' + figmid + '!' +           \
              '!'.join(tfig.split('!')[2:])

            graphvals = 'time:' + Arange + ':' + Brange + ':' + timekind + ';' +     \
              timefmt + ':' + timeaxis + ':' + ','.join(CFvplot) + ':' + Aline + ',' \
              + Bline + ':2.,2.:' + Akind + ',' + Bkind + ':2.,2.:,:' + figtit + ':' \
              + timelabel + ':0|auto:' +  finame.replace('.'+kfig, '') + ':' +       \
              kfig + ':True' 

            fvarS = ','.join(vplot)

            drwins = 'draw_2lines_time'
            plotins = 'python ' + pyH + '/drawing.py -f ' + figfs +' -o ' + drwins + \
              " -S '" + graphvals + "' -v " + fvarS

            try:
                with gen.Capturing() as output:
                    sout = sub.call(plotins, shell=True)
            except:
                print errmsg
                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
                print sout
                for s1out in output: print s1out
                quit(-1)

            # keeping all figures
            trkobjf.write('\n')
            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
            trkobjf.write(plotins + '\n')

        elif kplot == 'Nlines':
            figtit = tfig.replace('!','|')

            linestdn = CFvplot[0]
            figfs = ','.join(fplot)

            Nlines = len(fplot)

            lablines = []
            for il in range(Nlines):
                linev = pplot[il]
                # Values are changed from specific values in `model_graphics.dat'
                if linev[3].find('%') == -1:
                    line = linev[2]
                    kind = linev[3]
                    mark = linev[4]
                    size = linev[5]
                    if il == 0:
                        vrange = str(linev[0]) + ',' + str(linev[1])
                        kinds = kind
                        lines = line
                        marks = mark
                        sizes = str(size)
                    else:
                        kinds = kinds + ',' + kind
                        lines = lines + ',' + line
                        marks = marks + '@' + mark
                        sizes = sizes + ',' + str(size)
                else:
                    lines = 'None'
                    kinds = '-'
                    marks = ','
                    sizes = '2.'

                secsf = fplot[il].split('/')
                Nsecsf = len(secsf)
                expvs = expgraphc[secsf[Nsecsf-3] + '/' + secsf[Nsecsf-2]]
                lablines.append(expvs.label)

            # It is assumed that if the space variable is 'lon': is desired a 
            #   (lon, vals) plot if it is 'lat': then (vals, lat) plot
            if gen.searchInlist(dplot[0],'lon'):
                spacedim ='lon'
                sdlab ='lon ($degrees\ East$)'
                axisvals = 'y'
                figmid = 'longitudinal'
            elif gen.searchInlist(dplot[0],'lat'):
                spacedim = 'lat'
                sdlab = 'lat ($degrees\ North$)'
                axisvals = 'x'
                figmid = 'meridional'
            elif gen.searchInlist(dplot[0],'pres'):
                spacedim = 'pres'
                sdlab = 'pres ($Pa$)'
                axisvals = 'x'
                figmid = 'vertical'
            else:
                print errmsg
                print '  ' + fname + ": in 'Nlines' only ready for: 'lon', 'lat'," + \
                 " 'pres' as common dimension !!"
                print '  dimensions in plot:', dplot[0]
                quit(-1)

            figtit = '!'.join(tfig.split('!')[0:2]) + '!' + figmid + '!' +           \
              '!'.join(tfig.split('!')[2:])
            figtit = figtit.replace('!',' ')
            leglabels = ','.join(lablines)

            graphvals = spacedim + ':' + axisvals + ':' + sdlab + ':auto:' +         \
              leglabels + ':' + CFvplot[0] +':'+ figtit + ':0|auto:' + lines + ':' + \
              kinds + ':' + marks +':'+ sizes +':'+ sizes +':all:'+                  \
              finame.replace('.'+kfig, '') + ':' + kfig + ':True'

            fvarS = vplot[0]

            drwins = 'draw_lines'
            plotins = 'python ' + pyH + '/drawing.py -f ' + figfs +' -o ' + drwins + \
              " -S '" + graphvals + "' -v " + fvarS

            try:
                with gen.Capturing() as output:
                    sout = sub.call(plotins, shell=True)
            except:
                print errmsg
                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
                print sout
                for s1out in output: print s1out
                quit(-1)

            # keeping all figures
            trkobjf.write('\n')
            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
            trkobjf.write(plotins + '\n')

        elif kplot == 'Nlines_time':
            figtit = tfig.replace('!','|')

            linestdn = CFvplot[0]
            figfs = ','.join(fplot)

            Nlines = len(fplot)

            lablines = []
            for il in range(Nlines):
                linev = pplot[il]
                # Values are changed from specific values in `model_graphics.dat'
                if linev[3].find('%') == -1:
                    line = linev[2]
                    kind = linev[3]
                    mark = linev[4]
                    size = linev[5]
                    if il == 0:
                        vrange = str(linev[0]) + ',' + str(linev[1])
                        kinds = kind
                        lines = line
                        marks = mark
                        sizes = str(size)
                    else:
                        kinds = kinds + ',' + kind
                        lines = lines + ',' + line
                        marks = marks + '@' + mark
                        sizes = sizes + ',' + str(size)
                else:
                    lines = 'None'
                    kinds = '-'
                    marks = ','
                    sizes = '2.'

                secsf = fplot[il].split('/')
                Nsecsf = len(secsf)
                expvs = expgraphc[secsf[Nsecsf-3] + '/' + secsf[Nsecsf-2]]
                lablines.append(expvs.label)

            valsaxis = 'y'
            figmid = 'evolution'

            figtit = '!'.join(tfig.split('!')[0:2]) + '!' + figmid + '!' +           \
              '!'.join(tfig.split('!')[2:])
            figtit = figtit.replace('!',' ')
            leglabels = ','.join(lablines)

            graphvals = 'time;' + valsaxis + ';time;' + leglabels + ';' +            \
              CFvplot[0] + ';' + figtit + ';None;time|' + '|'.join(tvals) +          \
              ';0|auto;' + kfig + ';' + kinds + ';' + lines + ';' + marks + ';' +    \
              sizes + ';' + sizes + ';all;-1;True' 

            fvarS = vplot[0]

            drwins = 'draw_lines_time'
            plotins = 'python ' + pyH + '/drawing.py -f ' + figfs +' -o ' + drwins + \
              " -S '" + graphvals + "' -v " + fvarS

            try:
                with gen.Capturing() as output:
                    sout = sub.call(plotins, shell=True)
            except:
                print errmsg
                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
                print sout
                for s1out in output: print s1out
                quit(-1)

            sout = sub.call('mv lines_time.' + kfig + ' ' + finame, shell=True)

            # keeping all figures
            trkobjf.write('\n')
            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
            trkobjf.write(plotins + '\n')

        elif kplot == 'shadcont2Dsfc':
            figtit = tfig.replace('!','|')
            shdstdn = CFvplot[0]
            cntstdn = CFvplot[1]
            figfs = ','.join(fplot)

            shad = pplot[0]
            srange = str(shad[0]) + ',' + str(shad[1])
            cbar = shad[2]
            cnt = pplot[1]
            crange = str(cnt[0]) + ',' + str(cnt[1])
            cfmt = cnt[3]
            ckind = cnt[4]
            cline = cnt[5]

            graphvals = ','.join(CFvplot)
            graphvals = graphvals + ':lon|-1,lat|-1:lon|-1,lat|-1:lon:lat:auto:' +   \
              cbar + ',auto,auto:' + ckind + ',' + cline + ':' + cfmt + ':' +        \
              srange + ':' + crange + ',9:' + figtit + ':' + kfig + ':None:' +       \
              mapval + ':True'

            fvarS = ','.join(vplot)

            drwins = 'draw_2D_shad_cont'
            plotins = "python " + pyH + "/drawing.py -f " + figfs + " -o " + drwins +\
              " -S '" + graphvals + "' -v " + fvarS
            try:
                with gen.Capturing() as output:
                    sout = sub.call(plotins, shell=True)
            except:
                print errmsg
                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
                print sout
                for s1out in output: print s1out
                quit(-1)

            sout = sub.call('mv 2Dfields_shadow-contour.'+kfig+' '+finame, shell=True)

            # keeping all figures
            trkobjf.write('\n')
            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
            trkobjf.write(plotins + '\n')

        elif kplot == 'shadconthovmsfc':
            figtit = tfig.replace('!','|')
            shdstdn = CFvplot[0]
            cntstdn = CFvplot[1]
            figfs = ','.join(fplot)

            shad = pplot[0]
            srange = str(shad[0]) + ',' + str(shad[1])
            cbar = shad[2]
            cnt = pplot[1]
            crange = str(cnt[0]) + ',' + str(cnt[1])
            cfmt = cnt[3]
            ckind = cnt[4]
            cline = cnt[5]
            # It is assumed that if the space variable is 'lon': is desired a 
            #   (lon, time) plot if it is 'lat': then (time, lat) plot
            if gen.searchInlist(dplot,'lon'):
                spacedim ='lon'
                figmid = 'longitudinal|evolution|of'
                figtit = '|'.join(tfig.split('!')[0:2]) + '|' + figmid + '|' +       \
                  '|'.join(tfig.split('!')[2:])
                reverse= 'None'
                dims= spacedim+'|-1,time|-1;'+spacedim+'|-1,time|-1;'+spacedim+';time'
            else:
                spacedim='lat'
                figmid = 'meridional|evolution|of'
                figtit = '|'.join(tfig.split('!')[0:2]) + '|' + figmid + '|' +       \
                  '|'.join(tfig.split('!')[2:])
                reverse='None'
                dims= spacedim+'|-1,time|-1;'+spacedim+'|-1,time|-1;time;'+spacedim

            graphvals = ','.join(CFvplot)
            graphvals = graphvals + ';' + dims +';auto;' + cbar + ',auto,auto' + ';' \
              + ckind + ',' + cline + ';' + cfmt + ';' + srange + ';' + crange +     \
              ',9;' + figtit + ';' + kfig + ';' + reverse + ';' + 'time|' + tunits + \
              '|'+ timekind + '|' + timefmt + '|' + timelabel + ';True'

            fvarS = ','.join(vplot)

            drwins = 'draw_2D_shad_cont_time'
            plotins = 'python ' + pyH + '/drawing.py -f ' + figfs +' -o ' + drwins + \
              " -S '" + graphvals + "' -v " + fvarS

            try:
                with gen.Capturing() as output:
                    sout = sub.call(plotins, shell=True)
            except:
                print errmsg
                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
                print sout
                for s1out in output: print s1out
                quit(-1)

            sout = sub.call('mv 2Dfields_shadow-contour.'+kfig+' '+finame, shell=True)

            # keeping all figures
            trkobjf.write('\n')
            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
            trkobjf.write(plotins + '\n')

        elif kplot == 'shadcont2Dzsec':
            figtit = tfig.replace('!','|')
            shdstdn = CFvplot[0]
            cntstdn = CFvplot[1]
            figfs = ','.join(fplot)

            shad = pplot[0]
            srange = str(shad[0]) + ',' + str(shad[1])
            cbar = shad[2]
            cnt = pplot[1]
            crange = str(cnt[0]) + ',' + str(cnt[1])
            cfmt = cnt[3]
            ckind = cnt[4]
            cline = cnt[5]

            # It is assumed that if the space variable is 'lon': is desired a 
            #   (lon, pres) plot it it is 'lat': then (pres, lat) plot
            if gen.searchInlist(dplot, 'lon'):
                spacedim='lon'
                figmid = 'long.|vert.|cross|sec.|of'
                figtit = '|'.join(tfig.split('!')[0:2]) + '|' + figmid + '|' +       \
                  '|'.join(tfig.split('!')[2:])          
                dims= spacedim+'|-1,pres|-1:'+spacedim+'|-1,pres|-1:'+spacedim+':pres'
            else:
                spacedim='lat'
                figmid = 'mer.|vert.|cross|sec.|of'
                figtit = '|'.join(tfig.split('!')[0:2]) + '|' + figmid + '|' +       \
                  '|'.join(tfig.split('!')[2:])
# NOT WORKING?   dims= spacedim+'|-1,pres|-1:'+spacedim+'|-1,pres|-1:'+'pres:'+spacedim
                dims= spacedim+'|-1,pres|-1:'+spacedim+'|-1,pres|-1:'+spacedim+':pres'

            reverse = 'None'

            graphvals = ','.join(CFvplot)
            graphvals = graphvals + ':' + dims + ':auto:' + cbar + ',auto,auto:' +   \
              ckind + ',' + cline + ':' + cfmt + ':' + srange + ':' + crange + ',9:' \
              + figtit + ':' + kfig + ':' + reverse + ':None:True'
 
            fvarS = ','.join(vplot)

            drwins = 'draw_2D_shad_cont'
            plotins = 'python ' + pyH + '/drawing.py -f ' + figfs + ' -o ' + drwins + \
              " -S '" + graphvals + "' -v " + fvarS
            try:
                with gen.Capturing() as output:
                    sout0 = sub.call(plotins, shell=True)
            except:
                print errmsg
                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
                for s1out in output: print s1out
                quit(-1)

            sout = sub.call('mv 2Dfields_shadow-contour.'+kfig+' '+finame, shell=True)

            # keeping all figures
            trkobjf.write('\n')
            trkobjf.write("# '" + kplot + "'; " + tfig.replace('!',' ') + '\n')
            trkobjf.write(plotins + '\n')

        else:
            print errmsg
            print "  " + fname + ": plot kind '" + kplot + "' not ready !!"
            quit(-1)

    trkobjf.close()

    if db: print "    End of '" + fname + "' "

    return

def opexplained(Opn,Exp):
    """ Function to provide the final explanation for the title in the graph. Some operations have 
      values separated by `~' which need to be processed
      Opn= full operation string
      Exp= explanation from dictionary
    >>> opexplained('dimvarvalue~pres~85000','@')
    pres@85000
    >>> opexplained('xmean','xmean')
    xmean
    """
    fname = 'opexplained'

    if Opn.find('~') == -1:
        explanation = Exp
    else:
        Opnv = Opn.split('~')
        if Opnv[0] == 'dimvarvalue':
            explanation =  Opnv[1] + Exp + Opnv[2]
        else:
            print errormsg
            print '  ' + fname + ": unknow '~' separated operation '" + Opnv[0] +    \
              "' !!"
            quit(-1)
    return explanation

def optitle(op,opexp):
    """ Function to creat the operation section in the title of a figure ('!' for spaces)
      op= '+' separated list of operations
      opexp= dictionary with the `explanation'(text in title) to appear for each operation
    >>> optitle('pinterp+xmean+tmean',{'pinterp':'pinterp', 'xmean':'xmean', 'tmean':'tmean'})
    $_{[pinterp\!xmean\!tmean]}$
    >>> optitle('pinterp+dimvarvalue~pres~85000+tmean',{'dimvarvalue':'@', 'pinterp':'pinterp', 'xmean':'xmean', 'tmean':'tmean'})
    $_{[pinterp\!pres@85000\!tmean]}$
    """
    fname = 'optitle'

    opvals = op.split('+')
    for op1 in opvals:
        op1n = op1.split('~')[0]
        if not opexp.has_key(op1n):
            print errmsg
            print '  ' + fname + ": no explanation for operation '" + op1n + "' !!"
            print '    provided:', opexp.keys()
            quit(-1)
        # Final explanation
        op1en = opexplained(op1,opexp[op1n]).replace(' ','!')

        if op1n == opvals[0].split('~')[0]:
            titop = '$_{[' + op1en
            if len(opvals) == 1: titop = titop + ']}$'
        elif op1n == opvals[len(opvals)-1].split('~')[0]:
            titop = titop + '\!' + op1en + ']}$'
        else:
            titop = titop + '\!' + op1en

    return titop

def create_figure_title(mod,exp,varops,explop):
    """ Function to create the title of a figure ('!' for spaces)
      mod: model name
      exp: experiment name
      varops: list with [var]|[op] values from wich create the title
      explop: dictionary with the `explanation'(text in title) to appear for each operation
    >>> expops = {'pinterp':'pinterp', 'xmean':'xmean', 'tmean':'tmean', 'last':'last'}
    >>> create_figure_title('WRF','current',['ua|pinterp+xmean+tmean', 'va|pinterp+xmean+tmean'], expops)
    WRF!current!ua!&!va$_{[pinterp\!xmean\!tmean]}$
    >>> create_figure_title('WRF','current',['ua|pinterp+xmean+tmean', 'ua|pinterp+xmean+last'], expops)
    WRF!current!ua$_{[pinterp\!xmean\!tmean]}$!&!!$_{[pinterp\!xmean\!last]}$
    >>> create_figure_title('WRF','current',['ua|pinterp+xmean+tmean', 'va|pinterp+xmean+last'], expops)
    WRF!current!ua$_{[pinterp\!xmean\!tmean]}$!&!va$_{[pinterp\!xmean\!last]}$
    >>> create_figure_title('WRF','current',['uas|dimvarvalue~lon~23.23', 'vas|dimvarvalue~lon~23.23'], expops)
    WRF!current!uas!&!vas$_{[lon@23.23]}$
    >>> create_figure_title('WRF','current',['hus|pinterp+tmean+dimvarvalue~pres~85000.', 'ta|pinterp+tmean+dimvarvalue~pres~85000.'], expops)
    WRF!current!hus!&!ta$_{[pinterp\!tmean\!pres@85000.]}$
    """
    fname = 'create_figure_title'

    titfig = mod + '!' + exp + '!'

    varns = []
    opns = []
    for varop in varops:
        vn = varop.split('|')[0]
        op = varop.split('|')[1]
        varns.append(vn)
        opns.append(op)
        if varop == varops[0]:
            allvarequal = True
            allopequal = True
            varinfig = vn
            opinfig = op
        else:
            if vn != varinfig: allvarequal = False
            if op != opinfig: allopequal = False

    Nvars = len(varns)
    Nops = len(opns)

    if allvarequal:
        titfig = titfig + '!' + varns[0]
        if allopequal:
            opS = optitle(op,explop)
            titfig = titfig + opS
        else:
            for opn in opns:
                opS = optitle(opn,explop)
                if opn == opns[0]:
                    titfig = titfig + opS 
                elif opn == opns[Nops-1]:
                    titfig = titfig + '!&!' + opS
                else:
                    titfig = titfig + ',!' + opS
    else: 
        if allopequal:
            opS = optitle(op,explop)
            for vn in varns:
                if vn == varns[0]:
                    titfig = titfig + '!' + vn
                elif vn == varns[Nvars-1]:
                    titfig = titfig + '!&!' + vn
                else:
                    titfig = titfig + ',!' + vn
            titfig = titfig + opS
        else:
            for ivop in range(Nvars):
                vn = varns[ivop]
                op = opns[ivop]
                vnop = vn + optitle(op,explop)
                
                if vn == varns[0]:
                    titfig = titfig + '!' + vnop
                elif vn == varns[Nvars-1]:
                    titfig = titfig + '!&!' + vnop
                else:
                    titfig = titfig + ',!' + vnop

    return titfig.replace('!!','!')
#expops = {'dimvarvalue':'@', 'last':'last', 'pinterp':'pinterp', 'xmean':'xmean', 'tmean':'tmean'}

def draw_plots(config, plots, mod, exp, odir, allvarcomp, figscr, debug):
    """ Function to draw all plots
      config= Configuration of the experiment
      plots= dictionary with the plots
      mod= model of the plots
      exp= experiment of the plots
      odir= output experiment folder
      allvarcomp= dictionary with all the variables to compute and their information
      figscr= whether figures should be done from the scratch or not

    * Plot as 
      {[kplot]} = [varn1]|[op1]#[varn2]|[op2]#[...[varnN]|[opN]], ...
        [kplot] ___ 
          diffmap2Dsfc: 2D map of surface differences values of 1 variable
          diffmap2Dz: 2D map of 3D differences values of 1 variable
          map2Dsfc: 2D map of surface values of 1 variable
          map3D: 2D map of 3D values of 1 variable
          shadconthovmsfc: Hovmoeller diagrams of 2 variable at the surface in shadow and the other in contourn
          shadcont2Dsfc: 2D map of shadow (1st variable) and countour (2nd variable) [stvar1]#[stvar2]
          shadcont2Dzsec: 2D map of vertical section of 2 variables one in shadow and the other in contourn
        [varn] 
          variable
        [op]
          '+' separated list of operations
        in figures with more than 1 variable, use '#' to separate the [varn]|[op]
    """ 
    fname = 'draw_plots'

    os.chdir(odir)

    # Dictionary with the operations with surnames for the operated variable
    opersurnames = {}
    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
    for opsr in opsur.keys():
        opn = opsr.split('_')[1]
        vls = opsur[opsr].split(':')
        opersurnames[opn] = vls

    # time values
    # Units time for the plots
    rd = config['CFreftime']
    tunits = config['CFunitstime'] + '!since!' + rd[0:4] + '-' + rd[4:6] + '-' +  \
      rd[6:8] + '!' + rd[8:10] + ':' + rd[10:12] + ':' + rd[12:14]
    # time ticks kind
    tkind = config['timekind']
    # time ticks format
    tfmt = config['timefmt']
    # time axis label
    tlab = config['timelabel']
    timevals = [tunits, tkind, tfmt, tlab]

    # Dictionary of plot specificities
    specplotkeyn = 'specificvarplot'
    plotspecifics = {}
    if config.has_key(specplotkeyn):
        #   [minval]: minimum value
        #   [maxval]: minimum value
        #   [colorbar]: name of the colorbar (from matplotlib) to use
        #   [cntformat]: format of the contour labels
        #   [colorcnt]: color for the countor lines
        plotspecs = config[specplotkeyn].split(':')
        for pltspc in plotspecs:
            pltvls = pltspc.split('|')
            vn = pltvls[0]
            op = pltvls[1]
            fn = pltvls[2]
            plotspecifics[fn + '_' + vn + '_' + op] = pltvls[3:]
        if debug:
            print 'Specific values for plots _______'
            gen.printing_dictionary(plotspecifics)

    # Kind of figures
    kindfigure = config['kindfig']

    # Graphical labels and configuration
    modgc, expgc = graphical_characteristics(config, debug)
    modv = modgc[mod]
    modlab = modv.label
    expv = expgc[mod+'/'+exp]
    explab = expv.label

    # Map value
    mapvalue = config['mapval']

    # pythone scripts HOME
    pyHOME = config['pyHOME']

    # Title-text of operations
    opexplained = {}
    optits = config['titleoperations'].split(':')
    for optit in optits:
        opn = optit.split('|')[0]
        opt = optit.split('|')[1]
        opexplained[opn] = opt
    if debug:
        print 'Titles for operations  _______'
        gen.printing_dictionary(opexplained)

    for kplot in plots.keys():
        varsplt = plots[kplot]
        for varplt in varsplt:
            if debug:
                print "  printing '" + kplot + "' var ':" + varplt + "'..."
            varops = varplt.split('#')

            # CF variables in plot
            CFvarsplot = []
            # Files in plot
            filesplot = []
            # Variables in plot within the files
            varsplot = []
            # Dims in figure
            dimsplot = []
            # pictoric values in figure
            pictplot = []
            # Name of the figure
            figname = ''
            # Title of the figure
            titfigure = ''

            ivp = 0
            for varop in varops:
                vn = varop.split('|')[0]
                op = varop.split('|')[1]

                # CF variables in plot
                CFvarsplot.append(vn)
 
                vnopS = vn + '_' + op
                if not allvarcomp.has_key(vnopS):
                    print errmsg
                    print '  ' + fname + ": no file for variable-operation '" +     \
                      vnopS + "' !!"
                vopvals = allvarcomp[vnopS]
                headf = vopvals[0]
                globalP = vopvals[3]
                gP = pinterpS(globalP)

                filen = vn + '_' + headf + gP + '_' + op.replace('+','_') + '.nc'
                filesplot.append(filen)
                # Do we have processed the given variable?
                if not os.path.isfile(filen):
                    print warnmsg
                    print "  " + fname + ": there is no file for variable '" + varop \
                      + "' skiping it !!"
                    break

                # Name of the variable inside the file
                vnsur = varnoper(vn, op, opersurnames)
                varsplot.append(vnsur)

                # Dimensions in file
                try:
                    with gen.Capturing() as output:
                        dims = ncvar.idims(filen)
                except:
                    print errmsg
                    print 'ncvar.idims('+filen+')'
                    for s1out in output: print s1out
                    quit(-1)

                dimsplot.append(dims)

                # pictoric values for the figure
                Sfivaop = kplot + '_' + vn + '_' + op
                if plotspecifics.has_key(Sfivaop):
                    pictvals = plotspecifics[Sfivaop]
                else:
                    Vvals = gen.variables_values(vn)
                    pictvals = [Vvals[2], Vvals[3], Vvals[6], '%g', 'fixc', 'black']

                pictplot.append(pictvals)

                # Header of the name of the figure
                if ivp == 0:
                    figname = kplot +'_'+ vn + '_' + headf + '_' + op.replace('+','-')
                else:
                    figname = figname + '-'+vn+'_'+headf+'_'+op.replace('+','-')

                ivp = ivp + 1
            # End of variable-operation
            figname = figname + '.' + kindfigure

            # Removig double '..' (from `dimvarvalue')
            figname = figname.replace('..', '.')

            # Title of figure
            titfigure = create_figure_title(modlab, explab, varops, opexplained)

            draw_plot(kplot, CFvarsplot, filesplot, varsplot, dimsplot, pictplot,    \
              figname, titfigure, kindfigure, mapvalue, timevals, expgc, odir,       \
              pyHOME, figscr, debug)

        # End of variables-operations

    # End of kind of plots

    if debug: print "    End of '" + fname + "' "

    return

def get_differences_var(vc,kdiff,db):
    """ Function to provide the operations to make for each variable from kdiff ('DIFFOP_' or 'DIFFVAR_') dictionary
      vc= kdiff dictionary, dictionary with all the parameters started with kdiff
      kdiff= kind of differences 'op', 'var'
    """
    fname = 'get_differences_var'

    if kdiff == 'op':
        diffn = 'DIFFOP'
    elif kdiff == 'var':
        diffn = 'DIFFVAR'
    else:
         print errmsg
         print '  ' + fname + ": differences '" + kdiff + "' not ready !!"
         quit(-1)
    LD = len(diffn + '_')

    iop = 1
    # list of operations by difference
    doopers = {}
    # operations by variable
    operationsvar = {}
    # individual operations by variable (not repeated)
    indivoperationsvar = {}

    # Variables with a calculation whic requires vertical interpolation to all the 
    #   file (operations stating by [kdiff]'_pinterp')
    LVp = len(diffn + '_pinterp')
    varglobalp = []

    for oper in vc.keys():
        Loper = len(oper)
        opn = oper[LD:Loper+1]
        vns = vc[oper].split(':')
        ops1 = opn.split('+')
        doopers[opn] = ops1
        for vn in vns:
            if not operationsvar.has_key(vn):
                operationsvar[vn] = [opn]
                indivoperationsvar[vn] = ops1
            else:
                opers = operationsvar[vn]
                opers.append(opn)
                operationsvar[vn] = opers
                opers1 = indivoperationsvar[vn]
                for op1 in ops1:
                    if not gen.searchInlist(opers1, op1):
                        opers1.append(op1)
                indivoperationsvar[vn] = opers1

        if oper[0:LVp] == diffn + '_pinterp':
            for vn in vns:
                if not gen.searchInlist(varglobalp,vn): varglobalp.append(vn)

        iop = iop + 1

    if db:
        print '  operations to make _______'
        gen.printing_dictionary(doopers)
        print '  operations by variable _______'
        gen.printing_dictionary(operationsvar)
        print '  individual operations by variable _______'
        gen.printing_dictionary(indivoperationsvar)
        print '  variables with the vertical interpolation for all the file _______'
        print '    #', varglobalp

    if db: print "    End of '" + fname + "' "

    return doopers, operationsvar, indivoperationsvar, varglobalp

def diffvarop_listconstruct(config, mods, exps, reprj, odir, debug):
    """ Function to create the list of differences to compute and draw
      config= Configuration of the experiment
      mods= list with the models to use
      exps= list with the experiments to use
      reprj= list with whether re-rpojection needs to be done
      odir= output differences folder
    """
    fname='diffvarop_listconstruct'
  
    diffop = gen.get_specdictionary_HMT(config, H='DIFFOP_')
    diffvar = gen.get_specdictionary_HMT(config, H='DIFFVAR_')

    if debug:
        if diffop is not None:
            print "  'op' differences ________"
            gen.printing_dictionary(diffop)
        if diffvar is not None:
            print "  'var' differences ________"
            gen.printing_dictionary(diffvar)

    # Getting variable characteristics from each model/experiment
    idir1 = config['ofold'] + '/' + mods[0] + '/' + exps[0]
    idir2 = config['ofold'] + '/' + mods[1] + '/' + exps[1]

    fvarcomp = 'varcompute.inf'
    if reprj[0]: fvc1 = 'reproj_' + fvarcomp
    else: fvc1 = fvarcomp

    if reprj[1]: fvc2 = 'reproj_' + fvarcomp
    else: fvc2 = fvarcomp

    files1, testfiles1, allcompvar1, Nvar1= read_varcomp_file(idir1 + '/' + fvc1)
    files2, testfiles2, allcompvar2, Nvar2= read_varcomp_file(idir2 + '/' + fvc2)

    diffops = {}
    if diffop is not None:
        # Getting 'op' differences characteristics
        doops, opvar, indivopvar, varglobalp = get_differences_var(diffop,'op',debug)

        Ndiffop = 0
        for vn in opvar.keys():
            vnops = opvar[vn]
            for op in vnops:
                vnop = vn + '_' + op
                if not allcompvar1.has_key(vnop):
                    print warnmsg
                    print '  ' + fname + ": DIFFOP variable+operation '" + vnop +    \
                      "' in '" + mods[0] + '/' + exps[0] + "' was not computed !!"
                    print '    skipping it!'
                    print '    variables computed:', allcompvar1
                    break
                if not allcompvar2.has_key(vnop):
                    print warnmsg
                    print '  ' + fname + ": DIFFOP variable+operation '" + vnop +    \
                      "' in '" + mods[1] + '/' + exps[1] + "' was not computed !!"
                    print '    skipping it!'
                    print '    variables computed:', allcompvar2
                    break
                vals1 = allcompvar1[vnop]
                vals2 = allcompvar2[vnop]

                headerf1 = vals1[0]
                headerf2 = vals2[0]

                diffops[vnop] = [mods[0],mods[1],exps[0],exps[1],headerf1,headerf2]

                Ndiffop = Ndiffop + 1        
    else:
        Ndiffop = 0

    diffvars = {}
    if diffvar is not None:
        # Getting 'var' differences characteristics
        doops,opvar,indivopvar,varglobalp = get_differences_var(diffvar,'var',debug)

        Ndiffvar = 0
        for vn in opvar.keys():
            vnops = opvar[vn]
            for op in vnops:
                vnop = vn + '_' + op
                if not allcompvar1.has_key(vnop):
                    print warnmsg
                    print '  ' + fname + ": DIFFVAR variable+operation '" + vnop +   \
                      "' in '" + mods[0] + '/' + exps[0] + "' was not computed !!"
                    print '    skipping it!'
                    print '    variables computed:', allcompvar1
                    break
                if not allcompvar2.has_key(vnop):
                    print warnmsg
                    print '  ' + fname + ": DIFFVAR variable+operation '" + vnop +   \
                      "' in '" + mods[1] + '/' + exps[1] + "' was not computed !!"
                    print '    skipping it!'
                    print '    variables computed:', allcompvar2
                    break
                vals1 = allcompvar1[vnop]
                vals2 = allcompvar2[vnop]
    
                headerf1 = vals1[0]
                headerf2 = vals2[0]

                diffvars[vnop] = [mods[0],mods[1],exps[0],exps[1],headerf1,headerf2]

                Ndiffvar = Ndiffvar + 1
    else:
        Ndiffvar = 0

    if diffop is not None:
        Sopdiffs = gen.dictKeysVals_stringList(diffops)
    else:
        Sopdiffs = 'none'
    # Outwritting the op diffferences file to avoid next time
    objf = open(odir + '/diffop.inf', 'w')
    objf.write('differences: ' + Sopdiffs + '\n')
    objf.write('Ndiff: ' + str(Ndiffop) + '\n')
    objf.close() 

    if diffvar is not None:
        Svardiffs = gen.dictKeysVals_stringList(diffvars)
    else:
        Svardiffs = 'none'
    # Outwritting the var diffferences file to avoid next time
    objf = open(odir + '/diffvar.inf', 'w')
    objf.write('differences: ' + Svardiffs + '\n')
    objf.write('Ndiff: ' + str(Ndiffvar) + '\n')
    objf.close() 

    if debug: print "    End of '" + fname + "' "

    return diffops, diffvars, Ndiffop, Ndiffvar

def read_diff_file(difffile):
    """ Function to read the file with the information about the differences
      difffile= file with the information
    """
    fname = 'read_diff_file'

    if not os.path.isfile(difffile):
        print errormsg
        print '  ' + fname + ": differences file '" + difffile + "' does not exist !!"
        quit(-1)

    objf = open(difffile, 'r')
    for line in objf:
        if line[0:1] != '#' and len(line) > 1:
            values = line.replace('\n','').split(' ')
            if values[0] == 'differences:':
                if values[1] != 'none':
                    diffs = gen.stringList_dictKeysVals(values[1])
                else:
                    diffs = {} 
            elif values[0] == 'Ndiff:':
                Ndiffs = int(values[1])

    objf.close()

    return diffs, Ndiffs

def compute_op_diffs(config, diffallop, odir, diffscr, debug):
    """ Function to compute operation differences
      config= experiment configuration
      diffallop= dictonary with the differences to compute
      odir= output folder
      diffscr= wether it should be done from the scratch
    """
    fname = 'compute_op_diffs'

    # output folder
    ofold = config['ofold']

    # Home of the python scripts
    pyH = config['pyHOME']

    #  opsur = operation surnames
    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
    Opsurs = {}
    for opk in opsur:
        opn = opk.split('_')[1]
        vals = opsur[opk]
        Opsurs[opn] = vals.split(':')

    # Getting in working dir
    os.chdir(odir)

    # CF dimensions
    CFdims = ['lon', 'lat', 'pres', 'time']

    # Trakcing file
    trkobjf = open(odir + '/all_diffop.inf', 'a')

    for vnop in diffallop.keys():
        vn = vnop.split('_')[0]
        op = vnop.split('_')[1]
        opS = op.replace('+','_')

        # `p' header add related to 'pinterp'
        if opS.find('pinterp') != -1: SgP = 'p'
        else: SgP = ''

        # Variable name within the file
        vninF = varnoper(vn, op, Opsurs)

        vals = diffallop[vnop]
        mod1 = vals[0]
        mod2 = vals[1]
        exp1 = vals[2]
        exp2 = vals[3]
        headerf1 = vals[4]
        headerf2 = vals[5]
        modexpdiff = mod2 + '/' + exp2 + ' - ' + mod1 + '/' + exp1
        modexpdiffS = mod2 + '-' + exp2 + '_' + mod1 + '-' + exp1

        ifile1 = ofold + '/' + mod1 + '/' + exp1 + '/' + vn + '_' + headerf1 + SgP + \
          '_' + opS + '.nc'
        ifile2 = ofold + '/' + mod2 + '/' + exp2 + '/' + vn + '_' + headerf2 + SgP + \
          '_' + opS + '.nc'
        difffilen = odir + '/' + vn + '_diffop_' + modexpdiffS + '_' + opS + '.nc'

        if diffscr:
            sout = sub.call('rm ' + difffilen + ' >& /dev/null', shell=True)

        if not os.path.isfile(difffilen):
            if debug:
                print "    Computing operation difference '" + op + "' of '" + vn +  \
                  "' for: '" + modexpdiff
                print "    var in file: '" + vninF + "'"

            # Range of the dimensions
            for CFd in CFdims:
                if CFd == CFdims[0]:
                    CFdimS = CFd + '|' + CFd + '|-1'
                else:
                    CFdimS = CFdimS + ';' + CFd + '|' + CFd + '|-1'

            ncvar.ivattrs(ifile2,vninF)

            # Attributes of the variable
            try:
                with gen.Capturing() as output:
                    varattrs = ncvar.ivattrs(ifile2,vninF)
            except:
                print errmsg
                print 'ncvar.ivarattrs(' + ifile2 + ',' + vninF + ')'
                for s1out in output: print s1out
                quit(-1)

            if debug:
                for s1out in output: print s1out
            varvals = vninF + ',' + varattrs['long_name'][0] + ',' +                 \
              varattrs['units'][0]

            values = CFdimS + '@' + 'add|' + ifile2 + '|' + vninF + '%' + CFdimS +   \
              '@' + 'sub|' + ifile1 + '|' + vninF
            pyins = 'python ' + pyH + "/nc_var.py -o compute_opersvarsfiles -S '" +  \
              values + "' -v '" + varvals + "'"

            try:
                with gen.Capturing() as output:
                    ncvar.compute_opersvarsfiles(values,varvals)
            except:
                print errmsg
                print 'ncvar.compute_opersvarsfiles(' + values + ',' + varvals + ')'
                for s1out in output: print s1out
                quit(-1)

            sout = sub.call('mv opersvarsfiles_'+vninF+'.nc ' + difffilen, shell=True)

            # keeping all differences
            trkobjf.write('\n')
            trkobjf.write("#" + vn + ' ' + op + '\n')
            trkobjf.write( pyins + '\n')

    trkobjf.close()
    if debug: print "    End of '" + fname + "' "

    return

def compute_diff_statistics(Ecnf, filen, ofold, varn, Opers, scr, db):
    """ Function to compute statistics from var diff files
      Ecnf= dictionary with the configuration of the experiment
      filen= variable difference file to use
      ofold= directory to write the output files
      varn= name of the variable
      Opers= kind of operation: (as possible multiple consecutive combination of operations separated by '+'
        [calc1]+[calc2] will compute first [calc1] and then [calc2] 
        acc: temporal accumulated values
        diff: differences between models
        dimvarvalue~[dimvarn]~[value]: Slicing along the index at the nearest given value of a dimension-variable
        direct: no statistics
        last: last temporal value
        Lmean: latitudinal mean values
        Lsec: latitudinal section (latitudinal value must be given, [var]@[lat]) 
        lmean: longitudinal mean values
        lsec: longitudinal section (longitudinal value must be given, [var]@[lat])
        pinterp: pressure interpolation (to the given $plevels)
        tmean: temporal mean values
        tstd: temporal standard deviation values
        tturb: Taylor's turbulence decomposition value (x - <x>) for time
        tvar: temporal variance values
        xmean: x-axis mean values
        xvar: x-axis variance values
        ymean: y-axis mean values
        zsum: vertical aggregated values
      scr= should it be done from the scratch?
    """
    fname='compute_diff_statistics'

    # Full header of file
    Lfilen = len(filen)
    Hfile = filen[0:Lfilen-3]

    # CF dimension-variables name
    OpersS = '+'.join(Opers)
    if OpersS.find('pinterp') != -1 or filen.find('pinterp') != -1:
        varnCFs = ['lon', 'lat', 'pres', 'time']
    else:
        varnCFs = ['lon', 'lat', 'time']

    # Variable-dimension CF dictionary
    CFdimvardict = {'lon': 'lon', 'lat': 'lat', 'pres': 'pres', 'time': 'time'}

    # Python scripts HOME
    pyH = Ecnf['pyHOME']

    #  opsur = operation surnames
    opsur = gen.get_specdictionary_HMT(Ecnf, H='opsur_',M='',T='')
    Opsurs = {}
    for opk in opsur:
        opn = opk.split('_')[1]
        vals = opsur[opk]
        Opsurs[opn] = vals.split(':')

    # Getting in working dir
    os.chdir(ofold)

    # File to keep track of all the statistics of var differences
    otrackf = open(ofold + '/all_vardiffstatistics.inf', 'a')

    if db: 
        print "    computing statistics of diff variable: '", varn ,"' operation '"+ \
      OpersS + "' using file '" + filen + " ...'"

    for op in Opers:
        # File name from previous operations, and name of the variable within the 
        #   file (some operations change the name of the variable)
        if op == Opers[0]:
            Fopers = op
            prevfile = filen
            vninF = varn
        else:
            Fopers = Fopers + '_' + op
            prevfile = fileon
        fileon = Hfile + '_' + Fopers + '.nc'

        SvarnCFs = ',' + ','.join(varnCFs) 
        CFvarnp = vninF + SvarnCFs

        if scr:
            sout = sub.call('rm ' + fileon + ' >& /dev/null', shell=True)

        # Just produce file in case output file does not exist
        if not os.path.isfile(fileon):
            if db: 
                print '  ' + fname + ": creation of file '" + fileon + "' ..."
            dofile = True
        else:
            dofile = False
        
        # Variables to be kept in the final file
        varkeep = []

        if op == 'acc':
            # temporal accumulated values
            print "  " + fname + ": kind '" + op + "' not ready !!"
            wrongstats.append(varn + '_' + opers)
            break

        elif op[0:11] == 'dimvarvalue':
            # Slicing along the index at the nearest given value of a dimension-variable
            dimvarn = op.split('~')[1]
            dimvalue = op.split('~')[2]
            vals = dimvarn + ',' + dimvalue + ',' + dimvalue + ',0'
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.DataSetSection_multivars(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'DataSetSection_multivars('+vals+',', prevfile, ',' +      \
                      CFvarnp + ')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv DataSetSection_multivars.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o DataSetSection_multivars -S '" + vals +   \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + CFvarnp + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing the given variable-dimension
            varnCFs.remove(dimvarn)

        elif op == 'direct':
            # no statistics
            if dofile:
                sout = sub.call('mv ' + prevfile + ' ' + fileon, shell=True)
        elif op == 'last':
            # last temporal value
            vals='time,-9,0,0'
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.DataSetSection(vals,prevfile)
                except:
                    print errmsg
                    print 'DataSetSection('+vals+',', prevfile, ')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                hprevfile = prevfile[0:len(prevfile)-3]

                sout = sub.call('mv ' + hprevfile + '_time_B-9-E0-I0.nc ' + fileon, \
                  shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o DataSetSection -S '" + vals + "' -f " +  \
                  prevfile
                otrackf.write("\n")
                otrackf.write("# " + varn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

        elif op =='Lmean':
            # latitudinal mean values
            print "  " + fname + ": kind '" + op + "' not ready !!"
            wrongstats.append(varn + '_' + opers)
            break
        elif op =='Lsec':
            # latitudinal section (latitudinal value must be given, [var]@[lat]) 
            print "  " + fname + ": kind '" + op + "' not ready !!"
            wrongstats.append(varn + '_' + opers)
            break
        elif op =='lmean':
            # longitudinal mean values
            print "  " + fname + ": kind '" + op + "' not ready !!"
            wrongstats.append(varn + '_' + opers)
            break
        elif op =='lsec':
            # longitudinal section (longitudinal value must be given, [var]@[lon])
            print "  " + fname + ": kind '" + op + "' not ready !!"
            wrongstats.append(varn + '_' + opers)
            break
        elif op == 'tmean':
            # temporal mean values
            vals='time|-1,time,mean,' + ':'.join(varnCFs)
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@')
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_mean.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + varn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'time'
            varnCFs.remove('time')

            varkeep.append('timestats')

        elif op == 'tstd':
            # temporal standard deviation values
            vals='time|-1,time,std,' + ':'.join(varnCFs) + ':' + vdnz
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_std.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + varn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'time'
            varnCFs.remove('time')

            varkeep.append('timestats')

        elif op == 'tturb':
            # temporal turbulent values
            vals='time|-1,time,turb,' + ':'.join(varnCFs)
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@')
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_turb.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + varn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            varkeep.append('timestats')

        elif op == 'tvar':
            # temporal variance values (Taylor's turbulence)
            vals='time|-1,time,var,' + ':'.join(varnCFs) + ':' + vdnz
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@') + dnz+'@'+vdnz
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_var.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + varn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'time'
            varnCFs.remove('time')

            varkeep.append('timestats')

        elif op == 'xmean':
            # x-axis mean values
            vals='lon|-1,lon,mean,' + ':'.join(varnCFs)
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@')
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_mean.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + varn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'lon'
            varnCFs.remove('lon')

            varkeep.append('lonstats')

        elif op == 'xvar':
            # x-axis variance values
            vals='lon|-1,lon,var,' + ':'.join(varnCFs)
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@')
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_var.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + varn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'lon'
            varnCFs.remove('lon')

            varkeep.append('lonstats')

        elif op == 'ymean':
            # y-axis mean values
            vals='lat|-1,lat,mean,' + ':'.join(varnCFs)
            dims = gen.dictvar_listS(varnCFs,CFdimvardict,',','@')
            if dofile:
                try:
                    with gen.Capturing() as output:
                        ncvar.file_oper_alongdims(vals,prevfile,CFvarnp)
                except:
                    print errmsg
                    print 'file_oper_alongdims('+vals+',', prevfile, ','+CFvarnp+')'
                    for s1out in output: print s1out
                    quit(-1)

                if db:
                    for s1out in output: print s1out

                sout = sub.call('mv file_oper_alongdims_mean.nc '+fileon, shell=True)

                # Keeping the operations
                pyins=pyH + "/nc_var.py -o file_oper_alongdims -S '" + vals +        \
                  "' -f " + prevfile + " -v " + CFvarnp
                otrackf.write("\n")
                otrackf.write("# " + varn + " " + Fopers + "\n")
                otrackf.write(pyins + "\n")

            # removing dimension variable-dimension 'lat'
            varnCFs.remove('lat')

            varkeep.append('latstats')

        elif op == 'zsum':
            # vertical aggregated values
            print "  " + fname + ": kind '" + op + "' not ready !!"
            wrongstats.append(varn + '_' + opers)
            break
        else:
            print errmsg
            print '  ' + fname + ": operation '" + op + "' not ready !!"
            quit(-1)

        # End of kind of operation

        # Variable name in file (vninF) changed due to operation 
        #   but only if previous operation does not the same 'statistic'
        chvn = gen.dictionary_key_list(Opsurs, op)
        if chvn is not None:
            oldvninF = vninF
            vninF = vninF + chvn
            CFvarnp = CFvarnp.replace(oldvninF,vninF)
                
        if dofile:
            if len(varkeep) > 0:
                varkeepS = ',' + ','.join(varkeep)
            else:
                varkeepS = ''

            # Adding such CF dimension variables which might not be in the 
            #   post-operation file
            try:
                with gen.Capturing() as output:  
                    varinfile = ncvar.ivars(fileon)
            except:
                print errmsg
                print 'ivars('+fileon+')'
                for s1out in output: print s1out
                # Removing file in order to make sure that it will be redone
                print '  ' + fname + " removing file '" + fileon + "' ..."
                sub.call('rm ' + fileon + ' >& /dev/null', shell=True)
                quit(-1)

            for CFvn in varnCFs:
                if not gen.searchInlist(varinfile, CFvn):
                    if db:
                        print '  ' + fname + ": recupering CF variable '" + CFvn + "'"
                    try:
                        with gen.Capturing() as output:
                            ncvar.fvaradd(prevfile+','+CFvn, fileon)
                    except:
                        print errmsg
                        print 'fvaradd(' + prevfile + ',' + CFvn +', ' + fileon + ')'
                        for s1out in output: print s1out
                        # Removing file in order to make sure that it will be redone
                        print '  ' + fname + " removing file '" + fileon + "' ..."
                        sub.call('rm ' + fileon + ' >& /dev/null', shell=True)
                        quit(-1)

                    if db:
                        for s1out in output: print s1out

            totalvarkeeps = CFvarnp + ',' + ','.join(varnCFs) + varkeepS
            try:
                with gen.Capturing() as output:
                    oclean = ncvar.cleaning_varsfile(totalvarkeeps,fileon)
            except:
                print errmsg
                print 'cleaning_varsfile('+totalvarkeeps+','+fileon+')'
                for s1out in output: print s1out
                # Removing file in order to make sure that it will be redone
                print '  ' + fname + " removing file '" + fileon + "' ..."
                sub.call('rm ' + fileon + ' >& /dev/null', shell=True)
                quit(-1)

            if db:
                for s1out in output: print s1out

    # End of operations
    otrackf.close()

    if db: print "    End of '" + fname + "' "

    return

def compute_var_diffs(config, diffallvar, odir, diffscr, debug):
    """ Function to compute variable differences
      config= experiment configuration
      diffallvar= dictonary with the differences to compute
      odir= output folder
      diffscr= wether it should be done from the scratch
    """
    fname = 'compute_var_diffs'

    # output folder
    ofold = config['ofold']

    # Home of the python scripts
    pyH = config['pyHOME']

    #  opsur = operation surnames
    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
    Opsurs = {}
    for opk in opsur:
        opn = opk.split('_')[1]
        vals = opsur[opk]
        Opsurs[opn] = vals.split(':')

    # Getting in working dir
    os.chdir(odir)

    # CF dimensions
    CFdims = ['lon', 'lat', 'pres', 'time']

    # Trakcing file
    trkobjf = open(odir + '/all_diffvar.inf', 'a')

    for vnop in diffallvar.keys():
        vn = vnop.split('_')[0]
        op = vnop.split('_')[1]
        opSfinal = op.replace('+','_')

        lop = op.split('+')
        # `p' header add related to 'pinterp'
        if op.find('pinterp') != -1:
            print warnmsg
            print '  ' + fname + ': there is the vertical interpolation as operation'
            print '    variable differences will be computed from that files where '
            print '    vertical interpolation was just done'
            ipop = lop.index('pinterp')
            opS = '_' + '_'.join(lop[0:ipop+1])
            doops = lop[ipop+1:]
            SgP = 'p'
            vninF = varnoper(vn, '+'.join(lop[0:ipop+1]), Opsurs)
            print '    operations to perform:', doops
        else: 
            SgP = ''
            opS = ''
            doops = lop
            vninF = vn

        Sdoops = '+'.join(doops)

        # Variable name within the file

        vals = diffallvar[vnop]
        mod1 = vals[0]
        mod2 = vals[1]
        exp1 = vals[2]
        exp2 = vals[3]
        headerf1 = vals[4]
        headerf2 = vals[5]
        modexpdiff = mod2 + '/' + exp2 + ' - ' + mod1 + '/' + exp1
        modexpdiffS = mod2 + '-' + exp2 + '_' + mod1 + '-' + exp1

        ifile1 = ofold + '/' + mod1 + '/' + exp1 + '/' + vn + '_' + headerf1 + SgP + \
          opS + '.nc'
        ifile2 = ofold + '/' + mod2 + '/' + exp2 + '/' + vn + '_' + headerf2 + SgP + \
          opS + '.nc'
        difffilen = odir + '/' + vn + '_diffvar_' + modexpdiffS + opS + '.nc'
        difffinalfilen = odir + '/' + vn + '_diffvar_' + modexpdiffS + opS + '_' +   \
          Sdoops + '.nc'

        if diffscr:
            sout = sub.call('rm ' + difffilen + ' >& /dev/null', shell=True)
            sout = sub.call('rm ' + difffinalfilen + ' >& /dev/null', shell=True)

        if not os.path.isfile(difffinalfilen):
            if debug:
                print "    Computing variable difference of '" + vn + "' for: '" +   \
                  modexpdiff + "' and performing operation '" + Sdoops + "'"
                print "    var in file: '" + vninF + "'"

            # Range of the dimensions
            for CFd in CFdims:
                if CFd == CFdims[0]:
                    CFdimS = CFd + '|' + CFd + '|-1'
                else:
                    CFdimS = CFdimS + ';' + CFd + '|' + CFd + '|-1'

            # Attributes of the variable
            with gen.Capturing() as output:
                varattrs = ncvar.ivattrs(ifile2,vninF)
            if debug: 
                for s1out in output: print s1out
            varvals = vninF + ',' + varattrs['long_name'][0] + ',' +                 \
              varattrs['units'][0]

            values = CFdimS + '@' + 'add|' + ifile2 + '|' + vninF + '%' + CFdimS +   \
              '@' + 'sub|' + ifile1 + '|' + vninF
            pyins = 'python ' + pyH + "/nc_var.py -o compute_opersvarsfiles -S '" +  \
              values + "' -v '" + varvals + "'"

            try:
                with gen.Capturing() as output:
                    ncvar.compute_opersvarsfiles(values,varvals)
            except:
                print errmsg
                print 'ncvar.compute_opersvarsfiles(' + values + ', ' + varvals + ')'
                for s1out in output: print s1out
                quit(-1)

            sout = sub.call('mv opersvarsfiles_'+vninF+'.nc ' + difffilen, shell=True)

            # keeping all differences
            trkobjf.write('\n')
            trkobjf.write("#" + vn + ' ' + op + '\n')
            trkobjf.write( pyins + '\n')

            # Computing statisitcs
            compute_diff_statistics(config, difffilen, odir, vninF, doops, diffscr,  \
              debug)

    trkobjf.close()

    if debug: print "    End of '" + fname + "' "

    return


def draw_diff_plots(config, plots, odir, allvarcomp, kdiff, figscr, debug):
    """ Function to draw all plots
      config= Configuration of the experiment
      plots= dictionary with the plots
      odir= output experiment folder
      allvarcomp= dictionary with all the variables to compute and their information
      kdiff= kind of differences:
        'diffop': plots from operation differences
        'diffvar': plots from variable differences
      figscr= whether figures should be done from the scratch or not

    * Plot as 
      {[kplot]} = [varn1]|[op1]#[varn2]|[op2]#[...[varnN]|[opN]], ...
        [kplot] ___ 
          diffmap2Dsfc: 2D map of surface differences values of 1 variable
          diffmap2Dz: 2D map of 3D differences values of 1 variable
          map2Dsfc: 2D map of surface values of 1 variable
          map3D: 2D map of 3D values of 1 variable
          shadconthovmsfc: Hovmoeller diagrams of 2 variable at the surface in shadow and the other in contourn
          shadcont2Dsfc: 2D map of shadow (1st variable) and countour (2nd variable) [stvar1]#[stvar2]
          shadcont2Dzsec: 2D map of vertical section of 2 variables one in shadow and the other in contourn
        [varn] 
          variable
        [op]
          '+' separated list of operations
        in figures with more than 1 variable, use '#' to separate the [varn]|[op]
    """ 
    fname = 'draw_diff_plots'

    os.chdir(odir)

    # Dictionary with the operations with surnames for the operated variable
    opersurnames = {}
    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
    for opsr in opsur.keys():
        opn = opsr.split('_')[1]
        vls = opsur[opsr].split(':')
        opersurnames[opn] = vls

    # time values
    # Units time for the plots
    rd = config['CFreftime']
    tunits = config['CFunitstime'] + '!since!' + rd[0:4] + '-' + rd[4:6] + '-' +  \
      rd[6:8] + '!' + rd[8:10] + ':' + rd[10:12] + ':' + rd[12:14]
    # time ticks kind
    tkind = config['timekind']
    # time ticks format
    tfmt = config['timefmt']
    # time axis label
    tlab = config['timelabel']
    timevals = [tunits, tkind, tfmt, tlab]

    if kdiff == 'diffop':
         specplotkeyn = 'specificdiffopplot'
    elif kdiff == 'diffvar':
         specplotkeyn = 'specificdiffvarplot'
    else:
         print errmsg
         print '  ' + fname + ": differences kind '" + kdiff + "' not ready !!"
         quit(-1)

    # Dictionary of plot specificities
    plotspecifics = {}
    if config.has_key(specplotkeyn):
        #   [minval]: minimum value
        #   [maxval]: minimum value
        #   [colorbar]: name of the colorbar (from matplotlib) to use
        #   [cntformat]: format of the contour labels
        #   [colorcnt]: color for the countor lines
        plotspecs = config[specplotkeyn].split(':')
        for pltspc in plotspecs:
            pltvls = pltspc.split('|')
            vn = pltvls[0]
            op = pltvls[1]
            fn = pltvls[2]
            plotspecifics[fn + '_' + vn + '_' + op] = pltvls[3:]
        if debug:
            print 'Specific values for plots _______'
            gen.printing_dictionary(plotspecifics)

    # Kind of figures
    kindfigure = config['kindfig']

    # Graphical labels and configuration
    modgc, expgc = graphical_characteristics(config, debug)

    # Map value
    mapvalue = config['mapval']

    # pythone scripts HOME
    pyHOME = config['pyHOME']

    # Title-text of operations
    opexplained = {}
    optits = config['titleoperations'].split(':')
    for optit in optits:
        opn = optit.split('|')[0]
        opt = optit.split('|')[1]
        opexplained[opn] = opt
    if debug:
        print 'Titles for operations  _______'
        gen.printing_dictionary(opexplained)

    for kplot in plots.keys():
        varsplt = plots[kplot]
        for varplt in varsplt:
            if debug:
                print "  printing '" + kplot + "' var ':" + varplt + "'..."
            varops = varplt.split('#')

            # CF variables in plot
            CFvarsplot = []
            # Files in plot
            filesplot = []
            # Variables in plot within the files
            varsplot = []
            # Dims in figure
            dimsplot = []
            # pictoric values in figure
            pictplot = []
            # Name of the figure
            figname = ''
            # Title of the figure
            titfigure = ''

            ivp = 0
            for varop in varops:
                vn = varop.split('|')[0]
                op = varop.split('|')[1]

                # CF variables in plot
                CFvarsplot.append(vn)
 
                vnopS = vn + '_' + op
                if not allvarcomp.has_key(vnopS):
                    print errmsg
                    print '  ' + fname + ": no file for variable-operation '" +     \
                      vnopS + "' !!"
                    print '  available ones:', allvarcomp.keys()
                    quit(-1)
                vopvals = allvarcomp[vnopS]
                mod1 = vopvals[0]
                mod2 = vopvals[1]
                mod3 = vopvals[2]
                mod4 = vopvals[3]
                headf = vopvals[4]

                modv = modgc[mod1]
                expv = expgc[mod1+'/'+exp1]
                modlab1 = modv.label
                explab1 = expv.label
                modv = modgc[mod2]
                expv = expgc[mod2+'/'+exp2]
                modlab2 = modv.label
                explab2 = expv.label

                modexpdiff = explab2 + '-' + explab1
                modexpdiffS = mod2 + '-' + exp2 + '_' + mod1 + '-' + exp1

                difffilen = odir + '/' + vn + '_' + kdiff + '_' + modexpdiffS +      \
                  '_' + op.replace('+','_') + '.nc'

                filesplot.append(difffilen)
                # Do we have processed the given difference?
                if not os.path.isfile(difffilen):
                    print warnmsg
                    print "  " + fname + ": there is no file for '" + kdiff +        \
                      "' difference '" + varop + "' skiping it !!"
                    break

                # Name of the variable inside the file
                vnsur = varnoper(vn, op, opersurnames)
                varsplot.append(vnsur)

                # Dimensions in file
                try:
                    with gen.Capturing() as output:
                        dims = ncvar.idims(difffilen)
                except:
                    print errmsg
                    print 'ncvar.idims('+difffilen+')'
                    for s1out in output: print s1out
                    quit(-1)

                dimsplot.append(dims)

                # pictoric values for the figure
                Sfivaop = kplot + '_' + vn + '_' + op
                if plotspecifics.has_key(Sfivaop):
                    pictvals = plotspecifics[Sfivaop]
                else:
                    Vvals = gen.variables_values(vn)
                    pictvals = [Vvals[2], Vvals[3], Vvals[6], '%g', 'fixsigc','black']

                pictplot.append(pictvals)

                # Header of the name of the figure
                if ivp == 0:
                    figname = kplot + '_' + vn + '_' + kdiff + '_' + modexpdiffS +   \
                      '_' + op.replace('+','-')
                else:
                    figname = figname + '-' + vn + '_' + op.replace('+','-')

                # Title of the figure:
                if mod1 == mod2:
                    modexptit = mod1 + '!' + exp2 + '-' + exp1
                else:
                    modexptit = modexpdiff

                ivp = ivp + 1
            # End of variable-operation
            figname = figname + '.' + kindfigure

            # Title of figure
            titfigure = create_figure_title(modexptit, kdiff, varops, opexplained)

            if len(titfigure) > 80:
                print warnmsg
                print ' ' + fname + ": figure title '" + titfigure.replace('!', ' ')+\
                  "' larger than 80 characters (actually", len(titfigure), ') !!'
                print "  simplifying it to '" + modexptit.replace('!',' ') + '!' +   \
                  kdiff + '!' + vnS + "'"
                titfigure = modexptit.replace(' ','!') + '!' + kdiff + '!' + vnS

            draw_plot(kplot, CFvarsplot, filesplot, varsplot, dimsplot, pictplot,    \
              figname, titfigure, kindfigure, mapvalue, timevals, expgc, odir,       \
              pyHOME, figscr, debug)

        # End of variables-operations

    # End of kind of plots

    if debug: print "    End of '" + fname + "' "

    return

def reproject_modexp(mod,exp,config,REPRJops,fscr,fadd,debug):
    """ Function to re-project and re-compute files from a given model/experiment
      mod= model
      exp= experiment
      config= configuration
      REPRJops= dictionary with the reprojection operators to use for reproject and for which variable
        #    remapbil   CDO: Bilinear interpolation
        #    remapbic   CDO: Bicubic interpolation
        #    remapdis   CDO: Distance-weighted average remapping
        #    remapnn    CDO: Nearest neighbor remapping
        # Operators only available if projections have the corners of the grid points.
        #    remapcon   CDO: First order conservative remapping
        #    remapcon2  CDO: Second order conservative remapping
        #    remaplaf   CDO: Largest area fraction remapping
        ####### ####### ####### or additionally ####### ####### #######
        # python remapping
        #    dis        python: Distance-weighted average remapping
        #    pnn        python: Nearest neighbor remapping

      fscr= re-projected files to compute from the scratch
      fadd= re-projected files to be added
    """
    fname = 'reproject_modexp'

    odir = config['ofold'] + '/' + mod + '/' + exp

    # NON re-projectable operations
    opNOreproj = config['NOreprojops'].split(':')

    # Need to pass to analyze all the data?
    inffiles = ['reproj_varcompute.inf', 'all_reproj_computevars.inf']
    if fscr:
        for inff in inffiles:
            if debug: print "    removing information file '" + inff + "' ..."
            ins = 'rm ' + odir + '/' + inff + ' >& /dev/null'
            sout = sub.call(ins, shell=True)

        objf = open(odir+'/all_reproj_computevars.inf','w')
        objf.write("## Computation of reprojected variables \n")
        objf.close()

    varcompf = odir + '/reproj_varcompute.inf'
    if fadd:
        sub.call('rm ' + varcompf +' >& /dev/null', shell=True)

    # CF dims-vardims pairs:
    CFvardims = ['long@lon', 'lat@lat', 'pres@pres', 'time@time']

    # Dictionary with the operations with surnames for the operated variable
    opersurnames = {}
    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
    for opsr in opsur.keys():
        opn = opsr.split('_')[1]
        vls = opsur[opsr].split(':')
        opersurnames[opn] = vls

    if not os.path.isfile(varcompf):
        # Getting original computing information
        origvarcompf = odir + '/varcompute.inf'
        origfiles,origtestfiles,origallcompvar, Nvar = read_varcomp_file(origvarcompf)
             
        Svarcompute = ''
        allcompvar = dict(origallcompvar)
        ivop = 0
        for vnop in origallcompvar.keys():
            vn = vnop.split('_')[0]
            op = vnop.split('_')[1]
            values = allcompvar[vnop]
            values[0] = 'reproj-' + values[0]
            allcompvar[vnop] = values

            # Building strings from list values
            model = values[1]
            diag = values[2]
            if model is not None:
                if type(model) == type(list([1, 2])):
                    Smodel = ':'.join(model)
                elif type(model) == type('A'):
                    Smodel = model
            else:
                Smodel = 'None'
            if diag is not None:
                if type(diag) == type(list([1, 2])):
                    Sdiag = ':'.join(diag)
                elif type(diag) == type('A'):
                    Sdiag = diag
            else:
                Sdiag = 'None'

            Svc = vn + '|' + op + '|' + values[0] + '|' + Smodel + '|' + Sdiag +     \
              '|' + values[3]
            if ivop == 0:
                Svarcompute = Svc 
            else:
                Svarcompute = Svarcompute + ',' + Svc 
            ivop = ivop + 1

        origstrfiles = gen.dictKeysVals_stringList(origfiles)
        origstrtestfiles = gen.dictKeysVals_stringList(origtestfiles)

        # Outwritting the varcompute to avoid next time (if it is not filescratch!)
        objf = open(odir + '/reproj_varcompute.inf', 'w')
        objf.write('files: ' + origstrfiles + '\n')
        objf.write('testfiles: ' + origstrtestfiles + '\n')
        objf.write('varcompute: ' + Svarcompute + '\n')
        objf.write('itotv: ' + str(ivop) + '\n')
        objf.close() 
    else: 
        if dbg:
            print warnmsg
            print '  ' + main + ": getting variables to compute already from file !!"
        origfiles, origtestfiles, allcompvar, Nvar = read_varcomp_file(varcompf)

    # End of avoiding to repeat all the experiment search

    print "    For experiment '" + exp + "' is required to re-project:", Nvar,       \
      "variables"

    if dbg:
        print 'Variables to re-project _______'
        gen.printing_dictionary(allcompvar)

    # Re-projection variables
    print "    Reprojecting variables ..."
    os.chdir(odir)

    # Keeping track of reprojection
    trckobjf = open(odir + '/all_reproj_computevars.inf','a')

    for vnop in allcompvar.keys():
        vn = vnop.split('_')[0]
        op = vnop.split('_')[1]
        values = allcompvar[vnop]
        VnOpS = varnoper(vn, op, opersurnames)

        # Only are reprojected that files which do contain 'lon' and 'lat', 
        #   otherways, it is only re-projected the original file
        inNO, iopNO, iop = gen.list_coincidences(opNOreproj,op.split('+'))

        reprjop = gen.dictionary_key_list(REPRJops,vn).split('_')[1]
        if reprjop is None:
            print errmsg
            print '  ' + fname + ": no reprojecting operator found for " +           \
              " variable '" + vn + "' !!"
            print '    available values _______'
            gen.printing_dictionary(REPRJops)
            quit(-1)
        
        # Re-projecting also the 'original' (at least with 'pinterp') file for varsdif
        if op.find('pinterp') != -1:
            ifileorig = vn +'_'+ values[0].replace('reproj-','') + 'p_pinterp.nc'
            ofileorig = vn +'_'+ values[0] + 'p_pinterp.nc'
            ofilen = vn +'_'+ values[0] + 'p_pinterp.nc'
            ifileop = vn +'_'+ values[0].replace('reproj-','') + 'p_' +              \
              op.replace('+','_') + '.nc'
        else:
            ifileorig = vn +'_'+ values[0].replace('reproj-','') + '.nc'
            ofileorig = vn +'_'+ values[0] + '.nc'
            ifileop = vn +'_'+ values[0].replace('reproj-','') + '_' +               \
              op.replace('+','_') + '.nc'

        ofileop = vn + '_' + values[0] + '_' + op.replace('+','_') + '.nc'

        if fscr:
            sout = sub.call('rm ' + odir+'/'+ofileorig + ' >& /dev/null', shell=True)
            sout = sub.call('rm ' + odir+'/'+ofileop + ' >& /dev/null', shell=True)

        # List of files to re-project (first without operation)
        ifiles = [odir+'/'+ifileorig, odir+'/'+ifileop]
        # List of output files after re-projection (first without operation)
        ofiles = [odir+'/'+ofileorig, odir+'/'+ofileop]

        for iif in range(2):
            ifile = ifiles[iif]
            ofile = ofiles[iif]

            # Only projecting original file-variable if there are not 'lon' and 'lat'
            doreproj = True
            if iif == 1 and len(inNO) > 0: 
                doreproj = False
                if dbg:
                    print warnmsg
                    print '  ' + fname + ": operation '" + op + "' can not be " +    \
                      "reprojected due to the lack of 'lon' or 'lat' skipping it !!"
                    print "    NOTE: difference can only be computed as:"
                    print "      - 1st: compute variable's differences between models"
                    print "      - 2nd: Perform the operations"

                # Building up resources to compute re-projected statistics
                if dbg:
                    print '  '+fname+': computing from the original reprojected file'
                ModI = ModelInf(mod, mod, 'lon', 'lat', 'pres', 'time', 'depth',     \
                  'lon', 'lat', 'pres', 'time', 'depth', None)
                if type(values) == type([1]): 
                  fileh = values[0]
                else:
                  fileh = values

                if type(values) == type([1]): 
                  fileh = values[0]
                else:
                  fileh = values

                VarI = VariableInf(vn, values[0], values[2], values[3])
                usefs = {values[0]: ofileorig}

                compute_statistics(ModI, config, config['ifold'], usefs,             \
                  config['ofold']+'/'+mod+'/'+exp,VarI,pinterp_var(op), op, fscr, dbg)

            if not os.path.isfile(ofile) and doreproj:
                if debug and ofile == ofiles[0]:
                    print "      reprojecting '" + vn + "' '" + op + "' into '" +    \
                      ofile + "'"

                # Reprojecting...
                if reprjop[0:6] == 'remap':
                    ins = config['cdoHOME'] + '/cdo ' + reprjop + ',' +              \
                      config['RefProjfile'] + ' ' + ifile + ' ' + ofile
                    try:
                        with gen.Capturing() as output:
                            sout = sub.call(ins, shell=True)
                    except:
                        print errmsg
                        print '  ' + fname + ': $ ' + ins
                        for s1out in output: print s1out
                        print sout
                        quit(-1)
                    if sout != 0:
                        print errmsg
                        print '  ' + fname + ': $ ' + ins
                        sub.call(ins, shell=True)
                        quit(-1)
                else:
                    # Name of the variable to reproject
                    if iif == 0:
                        reprjvn = vn
                    else: 
                        reprjvn = VnOpS
              
                    reprojvalues = 'lon,lat,' + config['RefProjfile'] + ',lon,lat,'+ \
                      reprjop + ',' + ':'.join(CFvardims)
                    ins= 'ncvar.reproject('+reprojvalues + ', ' + ifile + ', '+      \
                      reprjvn + ')'
                    try:
                        with gen.Capturing() as output:
                            sout=ncvar.reproject(reprojvalues, ifile, reprjvn)
                    except:
                        print errmsg
                        print '  ' + fname + ': ' + ins
                        for s1out in output: print s1out
                        quit(-1)

                    sub.call('mv reproject.nc ' + ofile, shell=True)

                if debug:
                    for s1out in output: print s1out
                    print sout

                # Keeping track of re-projection
                trckobjf.write('\n')
                if iif == 0:
                    trckobjf.write('# ' + vn + '\n')
                else:
                    trckobjf.write('# ' + vn + ' ' + op + '\n')

                trckobjf.write(ins + '\n')

                if reprjop[0:6] == 'remap':
                    # CFing the reprojected file
                    try:
                        with gen.Capturing() as output:
                            ncvar.CDO_toCF(ofile)
                    except:
                        print errmsg
                        print 'ncvar.CDO_toCF(' + ofile + ')'
                        for s1out in output: print s1out
                        # Removing file in order to make sure that it will be redone
                        print '  ' + fname + " removing file '" + ofile + "' ..."
                        sub.call('rm ' + ofile + ' >& /dev/null', shell=True)
                        quit(-1)

                    if debug:
                        for s1out in output: print s1out
                
    trckobjf.close()

    if debug: print "    End of '" + fname + "' "

    return

def allmodexps_listconstruct(config, modexps, odir, debug):
    """ Function to create the list of all model-experiments to draw
      config= Configuration of the experiment
      modexps= list with the models-experiments pairs to use
      odir= output differences folder
    """
    fname='allmodexps_listconstruct'
  
    allplts = gen.get_specdictionary_HMT(config, H='PLOTALLMODEXP_')

    if debug:
        if allplts is not None:
            print "  'all mod-exps' plots ________"
            gen.printing_dictionary(allplts)

    # Getting Model-Reference Projection
    ModComProj = cnf['RefProj'].split(':')

    allmodexp = {}
    Nallmodexp = 0
    for modexp in modexps:
        if debug:
            print "    " +  fname + ": '" + modexp + "' ..."
        mod = modexp.split('/')[0]
        exp = modexp.split('/')[1]

        # Getting variable characteristics from each model/experiment
        idir = config['ofold'] + '/' + modexp

        reproj = not gen.searchInlist(ModComProj,mod)

        fvarcomp = 'varcompute.inf'
        if reproj: fvc = 'reproj_' + fvarcomp
        else: fvc = fvarcomp

        files, testfiles, allcompvar, Nvar= read_varcomp_file(idir + '/' + fvc)

        if allplts is not None:
            # Getting all model-experiments plots characteristics

            for allplt in allplts.keys():
                pltn = allplt.split('_')[1]
                if debug and modexp == modexps[0]:
                    print '    ' + fname + ":      plot '" + pltn + "' "
                varops = allplts[allplt].split(':')
                for varop in varops:
                    vn = varop.split('|')[0]
                    op = varop.split('|')[1]
                    vnop = vn + '_' + op
                    if debug and modexp == modexps[0]:
                        print '    ' + fname + ":        varop '" + vnop + "' "
                    if not allcompvar.has_key(vnop):
                        print warnmsg
                        print '  ' + fname + ": variable+operation '" + vnop +       \
                          "' in '" + modexp + "' was not computed !!"
                        print '    skipping it!'
                        print '    computed values:', allcompvar.keys()
                        break

                    vals = allcompvar[vnop]

                    headerf = vals[0]
                    if modexp == modexps[0]:
                        allmodexp[vnop] = [modexp,headerf]
                        Nallmodexp = Nallmodexp + 1
                    else:
                        prevals = allmodexp[vnop]
                        prevals = prevals + [modexp,headerf]
                        allmodexp[vnop] = prevals

        else:
            Nallmodexp = 0

    if allplts is not None:
        Sallmodexps = gen.dictKeysVals_stringList(allmodexp)
    else:
        Sallmodexps = 'none'

    # Outwritting the all model-experiment plots file to avoid next time
    objf = open(odir + '/allmodexp.inf', 'w')
    objf.write('allmodexps: ' + Sallmodexps + '\n')
    objf.write('Nallmodexp: ' + str(Nallmodexp) + '\n')
    objf.close()

    if debug: print "    End of '" + fname + "' "

    return allmodexp, Nallmodexp

def allmodexps_read(filen):
    """ Function to read the file with the information about the all model-experiment plots
      filen= file with the information
    """
    fname = 'allmodexps_read'

    if not os.path.isfile(filen):
        print errormsg
        print '  ' + fname + ": all model-experiment file '" + filen +               \
          "' does not exist !!"
        quit(-1)

    objf = open(filen, 'r')
    for line in objf:
        if line[0:1] != '#' and len(line) > 1:
            values = line.replace('\n','').split(' ')
            if values[0] == 'allmodexps:':
                allplts = gen.stringList_dictKeysVals(values[1])
            elif values[0] == 'Nallmodexp:':
                Nallmodexp = int(values[1])

    objf.close()

    return allplts, Nallmodexp

def draw_allmodexp_plots(config, allvarcomp, plots, modexp, odir, figscr, debug):
    """ Function to draw all model-experiment plots
      config= Configuration of the experiment
      allvarcomp = dictionary with the file headers for each variable
      plots= dictionary with the plots
      modexps= list with the names of the model-experiment pairs
      odir= output experiment folder
      figscr= whether figures should be done from the scratch or not

    * Plot as 
      {[kplot]} = [varn1]|[op1]#[varn2]|[op2]#[...[varnN]|[opN]], ...
        [kplot] ___ 
          Nlines: Multiple lines with one line by model-experiment pairs
        [varn] 
          variable
        [op]
          '+' separated list of operations
        in figures with more than 1 variable, use '#' to separate the [varn]|[op]
    """ 
    fname = 'draw_allmodexp_plots'

    os.chdir(odir)

    # Folder with the processed files
    infold = config['ofold']

    # Dictionary with the operations with surnames for the operated variable
    opersurnames = {}
    opsur = gen.get_specdictionary_HMT(config, H='opsur_',M='',T='')
    for opsr in opsur.keys():
        opn = opsr.split('_')[1]
        vls = opsur[opsr].split(':')
        opersurnames[opn] = vls

    # time values
    # Units time for the plots
    rd = config['CFreftime']
    tunits = config['CFunitstime'] + '!since!' + rd[0:4] + '-' + rd[4:6] + '-' +  \
      rd[6:8] + '!' + rd[8:10] + ':' + rd[10:12] + ':' + rd[12:14]
    # time ticks kind
    tkind = config['timekind']
    # time ticks format
    tfmt = config['timefmt']
    # time axis label
    tlab = config['timelabel']
    timevals = [tunits, tkind, tfmt, tlab]

    # Dictionary of plot specificities
    specplotkeyn = 'specificvarplot'
    plotspecifics = {}
    if config.has_key(specplotkeyn):
        #   [minval]: minimum value
        #   [maxval]: minimum value
        #   [colorbar]: name of the colorbar (from matplotlib) to use
        #   [cntformat]: format of the contour labels
        #   [colorcnt]: color for the countor lines
        plotspecs = config[specplotkeyn].split(':')
        for pltspc in plotspecs:
            pltvls = pltspc.split('|')
            vn = pltvls[0]
            op = pltvls[1]
            fn = pltvls[2]
            plotspecifics[fn + '_' + vn + '_' + op] = pltvls[3:]
        if debug:
            print 'Specific values for plots _______'
            gen.printing_dictionary(plotspecifics)

    # Graphical labels and configuration
    modgc, expgc = graphical_characteristics(config, debug)

    # Kind of figures
    kindfigure = config['kindfig']

    # Map value
    mapvalue = config['mapval']

    # pythone scripts HOME
    pyHOME = config['pyHOME']

    # Title-text of operations
    opexplained = {}
    optits = config['titleoperations'].split(':')
    for optit in optits:
        opn = optit.split('|')[0]
        opt = optit.split('|')[1]
        opexplained[opn] = opt
    if debug:
        print 'Titles for operations  _______'
        gen.printing_dictionary(opexplained)

    for kplot in plots.keys():
        varsplt = plots[kplot]
        for varplt in varsplt:
            if debug:
                print "  printing '" + kplot + "' var: '" + varplt + "'..."
            varops = varplt.split('#')

            # CF variables in plot
            CFvarsplot = []
            # Files in plot
            filesplot = []
            # Variables in plot within the files
            varsplot = []
            # Dims in figure
            dimsplot = []
            # pictoric values in figure
            pictplot = []
            # Name of the figure
            figname = ''
            # Title of the figure
            titfigure = ''

            ivp = 0
            for varop in varops:
                vn = varop.split('|')[0]
                op = varop.split('|')[1]

                # CF variables in plot
                CFvarsplot.append(vn)
 
                if op.find('pinterp') != -1: gP = 'p'
                else: gP = ''

                vnopS = vn + '_' + op                    
                if not allvarcomp.has_key(vnopS):
                    print errmsg
                    print '  ' + fname + ": no file for variable-operation '" +      \
                      vnopS + "' !!"
                modexpvals = allvarcomp[vnopS]
                Nmodexps = len(modexpvals)

                for imodexp in range(0,Nmodexps,2):
                    modexp = modexpvals[imodexp]
                    headf = modexpvals[imodexp+1]
                    if debug:
                        print '  ' + fname + "    model/experiment: '" +  modexp + "'"

                    filen = infold + '/' + modexp + '/' + vn +'_'+ headf+gP + '_' +  \
                      op.replace('+','_') + '.nc'

                    filesplot.append(filen)
                    # Do we have processed the given variable?
                    if not os.path.isfile(filen):
                        print warnmsg
                        print "  " + fname + ": there is no file for variable '" +   \
                          varop + "' skiping it !!"
                        break

                    if imodexp == 0:
                        # Name of the variable inside the file
                        vnsur = varnoper(vn, op, opersurnames)
                        varsplot.append(vnsur)

                    # Dimensions in file
                    try:
                        with gen.Capturing() as output:
                            dims = ncvar.idims(filen)
                    except:
                        print errmsg
                        print 'ncvar.idims('+filen+')'
                        for s1out in output: print s1out
                        quit(-1)

                    if imodexp == 0:
                        # Dimensions in plot
                        dimsplot.append(dims)

                        # Header of the name of the figure
                        if ivp == 0:
                            figname = kplot+'_'+vn+'_'+headf+'_'+ op.replace('+','-')
                        else:
                            figname = figname+'-'+vn+'_'+headf+'_'+op.replace('+','-')

                    # pictoric values for the figure
                    Sfivaop = kplot + '_' + vn + '_' + op
                    if plotspecifics.has_key(Sfivaop):
                        pictvals = plotspecifics[Sfivaop]
                    else:
                        Vvals = gen.variables_values(vn)
                        pictvals = [Vvals[2],Vvals[3],Vvals[6], '%g', 'fixc', 'black']

                    pictplot.append(pictvals)

                ivp = ivp + 1
            # End of variable-operation
            figname = figname + '.' + kindfigure

            # Title of figure
            titfigure = 'all!model-experiments!' + vn + optitle(op,opexplained)

            draw_plot(kplot, CFvarsplot, filesplot, varsplot, dimsplot, pictplot,    \
              figname, titfigure, kindfigure, mapvalue, timevals, expgc, odir,       \
              pyHOME, figscr, debug)

        # End of variables-operations

    # End of kind of plots

    if debug: print "    End of '" + fname + "' "

    return

def figinf_figname(figname,config):
    """ Function to retrieve figure information from its name
      figname= name of the figure
      config= configuration of `model_graphics.py'
    >>> figinf_figname('2linesTime_wss_wrfout_xvar-ymean-tas_wrfout_xvar-ymean.pdf', cnf)
    (2, '2linesTime', 'wss-tas', 'xvar-ymean@xvar-ymean')
    >>> figinf_figname('Nlines_tas_wrfout_tturb-xmean-tmean.pdf', cnf)
    (1, 'Nlines', 'tas', 'tturb-xmean-tmean')
    >>> figinf_figname('2lines_wss_diffvar_WRF-micro1_WRF-current_tturb-xmean-last-tas_tturb-xmean-last.pdf',cnf)
    (2, '2lines', 'wss-tas', 'tturb-xmean-last@tturb-xmean-last')
    """
    fname = 'figinf_figname'
    plotsecs = figname.split('_')
    plotk = plotsecs[0]

    if gen.searchInlist(config['pairfigures'].split(':'), plotk):
        Nvars = 2
        varn1 = plotsecs[1]
        if plotsecs[2] == 'diffvar' or plotsecs[2] == 'diffop':
            statn1varn2 = plotsecs[5].split('-')
            statn2 = plotsecs[6].split('.')[0]
            Lstatn1varn2 = len(statn1varn2)
            statn1 = '-'.join(statn1varn2[0:Lstatn1varn2-1])
            varn2 = statn1varn2[Lstatn1varn2-1]
        else:
            statn1varn2 = plotsecs[3].split('-')
            statn2 = plotsecs[5].split('.')[0]
            Lstatn1varn2 = len(statn1varn2)
            statn1 = '-'.join(statn1varn2[0:Lstatn1varn2-1])
            varn2 = statn1varn2[Lstatn1varn2-1]
        return Nvars, plotk, varn1+'-'+varn2, statn1+'@'+statn2
    elif gen.searchInlist(config['Nsourcesfigures'].split(':'), plotk):
        Nvars = 1
        varn = plotsecs[1]
        statn = plotsecs[3].split('.')[0]
        return Nvars, plotk, varn, statn
    else:
        print errmsg
        print '  ' + fname + ": kind of figure '" + plotk + "' not ready !!"
        print '    ready ones _______'
        print "    'pair_figures': figures of pair of variables: ",                  \
          config['pairfigures'].split(':')
        print "    'N-sources_figures': figures of pair of variables: ",                  \
          config['Nsourcesfigures'].split(':')
        quit(-1)
    
# Files with information about the configuration of the script
inffiles = ['varcompute.inf', 'all_computevars.inf', 'all_statsvars.inf']

#######    #######
## MAIN
    #######
if not os.path.isfile('model_graphics.dat'):
    print errmsg
    print main + ": configuration file 'model_graphics.dat' does not exist!!"
    quit(-1)

# Getting configuration from external ASCII file 'model_graphics.dat'
cnf = gen.get_configuration('model_graphics.dat', False)

verify_configuration(cnf, gen.Str_Bool(cnf['debug']))

# scratches
scratch, filescratch, figscratch, diffscratch, figdiffscratch, moddiffscratch,       \
  figmoddiffscratch, figallmodexpscratch, addfiles, addfigures, adddiffs,            \
  adddifffigures, addmoddiffs, addmoddifffigures, addallmodexpfigures, dbg =         \
  scratches(cnf)

if dbg:
    print 'Experiment configuration scratches _______'
    print '  scratch:', scratch, 'filescratch:', filescratch,'figscratch:', figscratch
    print '  diffscratch:', diffscratch, 'figdiffscratch:', figdiffscratch
    print '  moddiffscratch:', moddiffscratch, 'figmoddiffscratch:', figmoddiffscratch
    print 'Experiment configuration adding _______'
    print '  addfiles:', addfiles, 'addfigures:', addfigures
    print '  adddiffs:', adddiffs, 'adddifffigures:', adddifffigures
    print '  addmoddiffs:', addmoddiffs, 'addmoddifffigures:', addmoddifffigures

# Getting models
mods = cnf['models'].split(':')

# Getting graphical characeristics
modGC, expGC = graphical_characteristics(cnf, dbg)

# Models loop
##

# dictionary with the experiments of each model
modexps = {}

for mod in mods:
    print mod
    if scratch:
        if cnf['ofold'] != cnf['ifold']:
            sout = sub.call('rm -rf '+cnf['ofold']+'/'+mod+' > /dev/null', shell=True)
        else:
            print warnmsg
            print '  ' + main + ": input '" + cnf['ifold'] + "' and output '" +      \
              cnf['ofold'] + "' are the same folder !!"
            print "    Keeping output folder although it has 'scratch' start"

    # Get experiments and headers of model
    exps, fheaders = exp_headers(mod,cnf)

    modexps[mod] = exps
    # Characteristics of the model
    Modinf = ncvar.model_characteristics(mod,'None','False')
    dnx = Modinf.dimxn
    dny = Modinf.dimyn
    dnz = Modinf.dimzn
    dnt = Modinf.dimtn
    dns = Modinf.dimsn
    vdnx = Modinf.vardxn
    vdny = Modinf.vardyn
    vdnz = Modinf.vardzn
    vdnt = Modinf.vardtn
    vdns = Modinf.vardsn

    modgraphv = modGC[mod]
    Modinf = ModelInf(mod, mod, dnx, dny, dnz, dnt, dns, vdnx, vdny, vdnz, vdnt,     \
      vdns, modgraphv.tmod)

    if dbg:
        print '  model characteristics _______'
        gen.printing_class(Modinf)

    moddims = dnx + ',' + dny + ',' + dnz + ',' + dnt
    modvdims = vdnx + ',' + vdny + ',' + vdnz + ',' + vdnt

# Experiments loop
##
    for exp in exps:
        print '  ' + exp + '...'

        # input folder
        iwdir = cnf['ifold'] + '/' + mod + '/' + exp

        # Does input folder exist?
        if not os.path.isdir(iwdir):
            if filescratch*addfiles:
                print errmsg
                print "  " + main + ": folder '" + iwdir + "' does not exist !!"
                quit(-1)
            else:
                print warnmsg
                print "  " + main + ": folder '" + iwdir + "' does not exist !!"
                print "    but start of files from scratch 'filescratch':",          \
                  filescratch, "and there is not the add of any file 'addfiles':",   \
                  addfiles
                print '    so, keep moving'

        owdir = cnf['ofold'] + '/' + mod + '/' + exp
        sout = sub.call('mkdir -p ' + owdir, shell=True)

        # Need to pass to analyze all the data?
        if filescratch:
            for inff in inffiles:
                if dbg: print "    removing information file '" + inff + "' ..."
                ins = 'rm ' + owdir + '/' + inff + ' >& /dev/null'
                sout = sub.call(ins, shell=True)

            objf = open(owdir+'/all_computevars.inf','w')
            objf.write("## Computation of variables \n")
            objf.close()
            objf = open(owdir+'/all_statsvars.inf','w')
            objf.write("## Computation of statistics \n")
            objf.close()

        varcompf = owdir + '/varcompute.inf'
        if addfiles:
            if dbg:
                print '  ' + main + ": adding variables to compute removing " +      \
                  "file '" + varcompf + "' ..."
            sub.call('rm ' + varcompf +' >& /dev/null', shell=True)

        if not os.path.isfile(varcompf):
            # Does input folder has header files?
            ih=1
            # Dictionary with the list of files for each headers
            files = {}
            # Dictionary with a file fromor each headers
            testfiles = {}

            for fh in fheaders:
                if filescratch:
                    ins = 'rm '+ owdir+'/*_' + fh + '*.nc >& /dev/null'
                    sout = sub.call(ins, shell=True)
                files1h = gen.files_folder(iwdir,fh)
                if len(files1h) < 1:
                    print errmsg
                    print '  ' + main + ": folder '" + iwdir + "' does not " +   \
                      "contain files '" + fh + "*' !!"
                    quit(-1)
                files[fh] = files1h
                testfiles[fh] = files1h[0]

            if dbg:
                print '  Dictionary of files _______'
                gen.printing_dictionary(files)

                print '  Dictionary of test-files _______'
                gen.printing_dictionary(testfiles)
         
            allcompvar, Nvar = compvars_listconstruct(cnf, Modinf, files, testfiles, \
              iwdir, owdir, dbg)
        else:
            if dbg:
                print warnmsg
                print '  '+main+": getting variables to compute already from file !!"
            files, testfiles, allcompvar, Nvar = read_varcomp_file(varcompf)

        # End of avoiding to repeat all the experiment search

        print "    For experiment '" + exp + "' is required to compute:", Nvar,      \
          "variables"

        if dbg:
            print 'Variables to compute _______'
            gen.printing_dictionary(allcompvar)

        # Computing variables
        print "    Computing variables ..."
        compute_vars(cnf, Modinf, iwdir, owdir, files, allcompvar, filescratch, dbg)

        # Figures of variables direct from files
        print "    Plotting direct figures ..."
        dirfigf = owdir + '/directplotsdraw.inf'
        if figscratch:
            sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)

            objf = open(owdir+'/all_figures.inf','w')
            objf.write("## Drawing of all figures \n")
            objf.close()

        if addfigures:
            if dbg:
                print '  ' + main + ": adding direct figures removing file '" +      \
                  dirfigf + "' ..."
            sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)

        if not os.path.isfile(dirfigf):
            listplots, Nplt = plots_listconstruct(cnf, 'DIRPLT', dirfigf, owdir, dbg)
        else:
            if dbg:
                print warnmsg
                print '  ' + main + ": getting plots to draw already from file !!"
            listplots, Nplt = read_plot_file(dirfigf)

        # End of avoiding to repeat all the plots search

        print "    For experiment '" + exp + "' is required to plot:", Nplt, "plots"

        if dbg:
            print 'Plots to draw _______'
            gen.printing_dictionary(listplots)

        draw_plots(cnf, listplots, mod, exp, owdir, allcompvar, figscratch, dbg)

    # end of experiments loop

### ## #
# Experiment differences
## # ## #

# Experiments loop
##
    print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
    print "  ** '" + mod + "': Inter experiments differences  "
    print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
    print '    experiments:', exps
    # There are two kind of differences:
    #   DIFFop: differences between operations of each given variable. 
    #     Differences are computed directly from the last stage of the operation
    #   DIFFvar: operations of differences between variables
    #     First are computed the differences from the initial variable file
    #     and then operations are made
    # NOTE: remember that: meanvar2 - meanvar1 = mean(var2 - var1)
    difffiles = ['diffop.inf', 'diffvar.inf']

    Nexps = len(exps)
    for iexp1 in range(0,Nexps-1):
        exp1 = exps[iexp1]
        for iexp2 in range(iexp1+1,Nexps):
            exp2 = exps[iexp2]
            Sexps = exp2 + '-' + exp1
            print '  ' + Sexps + '...'
            owdir = cnf['ofold'] + '/' + mod + '/' + exp2 + '-' + exp1
            sout = sub.call('mkdir -p ' + owdir, shell=True)

            # Removing files with list of differences if should be started from scratch or add differences
            for fdiff in difffiles:
                difff = owdir + '/' + fdiff
                if diffscratch:
                    sub.call('rm -rf' + owdir +' >& /dev/null', shell=True)
                    sub.call('rm ' + difff +' >& /dev/null', shell=True)
                    objf = open(owdir+'/all_' + fdiff,'w')
                    objf.write("## Computation and drawing of differences " +        \
                      "between '" + exp2 + "'-'" + exp1 + "'\n")
                    objf.close()
                    difff = owdir + '/all_vardiffstatistics.inf'
                    sub.call('rm ' + difff +' >& /dev/null', shell=True)
                    objf = open(owdir+'/all_' + fdiff,'w')
                    objf.write("## Computation of all differences statistics " +     \
                      "between '" + exp2 + "'-'" + exp1 + "'\n")
                    objf.close()

                if adddiffs:
                    if dbg:
                        print '  ' + main + ": adding differences removing file '" + \
                          difff + "' ..."
                    sub.call('rm ' + difff +' >& /dev/null', shell=True)

            # op and var differences
            diffvarcompf = owdir + '/' + difffiles[0]
            if not os.path.isfile(diffvarcompf):
                alldiffop, alldiffvar,Nopdiffs,Nvardiffs=diffvarop_listconstruct(cnf,\
                  [mod, mod], [exp1, exp2], [False, False], owdir, dbg)
            else:
                if dbg:
                    print warnmsg
                    print '  ' + main + ": getting 'op' and 'var' differences to " + \
                      "compute already from file !!"
                alldiffop, Nopdiffops = read_diff_file(diffvarcompf)
                alldiffvar, Nopdiffvars = read_diff_file(owdir+'/'+difffiles[1])
                print 'Nopdiffops:', Nopdiffops,'Nopdiffvars:', Nopdiffvars

            # End of avoiding to repeat all the experiment search
            print "    For experiments '"+exp2+"'-'"+exp1+"' is required to " +      \
              "compute:", Nvar, "differences"

            if dbg:
                print 'Differences to compute _______'
                gen.printing_dictionary(alldiffop)
                gen.printing_dictionary(alldiffvar)

            # Computing differences
            ##
            print "    Computing operation differences ..."
            compute_op_diffs(cnf, alldiffop, owdir, diffscratch, dbg)
            print "    Computing variable differences ..."
            compute_var_diffs(cnf, alldiffvar, owdir, diffscratch, dbg)

            # Plotting operation differences
            ##
            print "  " + main + ": Plotting operation differences' figures ..."
            dirfigf = owdir + '/diffopplotsdraw.inf'
            if figscratch:
                sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)

                objf = open(owdir+'/all_diffopfigures.inf','w')
                objf.write("## Drawing of all operation difference figures \n")
                objf.close()

            if adddifffigures:
                if dbg:
                    print '  ' + main + ": adding differences' figures removing " +  \
                      "file '" + dirfigf + "' ..."

                sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)

            if not os.path.isfile(dirfigf):
                listplots, Nplt = plots_listconstruct(cnf, 'PLOTDIFFOP', dirfigf,    \
                  owdir, dbg)
            else:
                if dbg:
                    print warnmsg
                    print '  ' + main + ": getting plots to draw already from file !!"
                listplots, Nplt = read_plot_file(dirfigf)

            # End of avoiding to repeat all the plots search

            print "  For experiment 'operation' differences '" + Sexps + "' is " +   \
              "required to plot:" , Nplt, "plots"

            if dbg:
                print '    Plots to draw _______'
                gen.printing_dictionary(listplots)

            draw_diff_plots(cnf, listplots, owdir, alldiffop, 'diffop',              \
              figdiffscratch, dbg)

            # Plotting variable differences
            ##
            print "  " + main + ": Plotting variable differences' figures ..."
            dirfigf = owdir + '/diffvarplotsdraw.inf'
            if figscratch:
                sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)

                objf = open(owdir+'/all_diffvarfigures.inf','w')
                objf.write("## Drawing of all variables difference figures \n")
                objf.close()

            if adddifffigures:
                if dbg:
                    print '  '+main+": adding differences' figures removing file '" +\
                      dirfigf + "' ..."
                sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)

            if not os.path.isfile(dirfigf):
                listplots, Nplt = plots_listconstruct(cnf, 'PLOTDIFFVAR', dirfigf,   \
                  owdir, dbg)
            else:
                if dbg:
                    print warnmsg
                    print '  ' + main + ": getting plots to draw already from file !!"
                listplots, Nplt = read_plot_file(dirfigf)

            # End of avoiding to repeat all the plots search

            print "  For experiment 'variables' differences '" + Sexps + "' is " +   \
              "required to plot:" , Nplt, "plots"
 
            if dbg:
                print '    Plots to draw _______'
                gen.printing_dictionary(listplots)

            draw_diff_plots(cnf, listplots, owdir, alldiffvar, 'diffvar',            \
              figdiffscratch, dbg)

# end of mods loop

### ## #
# Model differences
### ## #

# Models already at the reference projection
ModComProj = cnf['RefProj'].split(':')
if dbg:
    print 'Models to use as reference for the projection:', ModComProj

# Operators to use for the reprojection of each variable
REPRJvar0 = gen.get_specdictionary_HMT(cnf, H='reprojectvar_')
REPRJvar = {}
for Skey in REPRJvar0.keys():
    REPRJvar[Skey] = gen.str_list(REPRJvar0[Skey],':')

if dbg:
    print 'Remaping operators to use _______'
    gen.printing_dictionary(REPRJvar)

# Getting common base weights for grid point remapping with CDO
frefgridn = 'RefProj.nc'
for mod in mods:
    frefgriddes = cnf['ofold'] + '/' + mod + '/' + frefgridn
    if moddiffscratch: 
        sout = sub.call('rm ' + frefgriddes + ' >& /dev/null', shell=True)
    if gen.searchInlist(ModComProj,mod) and not os.path.isfile(frefgriddes):
        # Looking for a file to use (all files are CF compilant, but we need to 
        #   be sure that at least they have 'lon', 'lat' variables)
        exps = modexps[mod]
        exp = exps[0]
        reprjfdir = cnf['ofold'] + '/' + mod + '/' + exp 
        reprjfile = None
        for reprjop in REPRJvar.keys():
            reprjvars = REPRJvar[reprjop]
            for reprjvar in reprjvars:
                reprjfs= gen.files_folder_HMT(folder=reprjfdir,head=reprjvar+'_',tail='.nc')
                for reprjf in reprjfs:
                    # Getting direct files as [var]_[fhead].nc
                    if len(reprjf.split('_')) == 2:
                        reprjfile = reprjfdir + '/' + reprjf
                        reprjvn = reprjvar
                        break
        if reprjfile is None:
            print errmsg
            print '  ' + main + ": no proper file to get projection information " +  \
              "from '" + reprjfdir + "' has been found !!"
            quit(-1)

        print "  Creation of reference projection information file from '" + mod+ "'"
        if dbg:
            print "    using file: '" + reprjfile

        # Getting only 1 time-step for disk space issues
        try:
             with gen.Capturing() as output:
                 ncvar.DataSetSection('time,0,1,1', reprjfile)
        except:
            print errmsg
            print 'ncvar.DataSetSection(time,0,1,1, ' + reprjfile + ')'
            for s1out in output: print s1out
            quit(-1)

        of = reprjfile.split('.')[0] + '_time_B0-E1-I1.nc'

        # Cleaning variables
        try:
             with gen.Capturing() as output:
                 ncvar.selvar('lon@lon,lat@lat,time@time', of, reprjvn)
        except:
            print errmsg
            print 'ncvar.selvar(lon@lon,lat@lat, ' + of + ', ' + reprjvn + ')'
            for s1out in output: print s1out
            quit(-1)
        sout = sub.call('mv selvar_new.nc ' + frefgriddes, shell=True)
        sout = sub.call('rm ' + of + ' >& /dev/null', shell=True)

        cnf['RefProjfile'] = frefgriddes

        if dbg:
            for s1out in output: print s1out
        break

    elif gen.searchInlist(ModComProj,mod) and os.path.isfile(frefgriddes):
        # Including reprojecting reference file in configuration dictionary
        cnf['RefProjfile'] = frefgriddes

Nmods = len(mods)
for imod1 in range(0,Nmods-1):
    mod1 = mods[imod1]
    exps1 = modexps[mod1]

    for imod2 in range(imod1+1,Nmods):
        mod2 = mods[imod2]
        exps2 = modexps[mod2]

        Smods = mod2 + '-' + mod1
        print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
        print "  ** '" + Smods + "': Inter models differences  "
        print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
        print '    experiments mod1:', exps1
        print '    experiments mod2:', exps2

# Experiments loop
##
        difffiles = ['diffop.inf', 'diffvar.inf']

        Nexps1 = len(exps1)
        Nexps2 = len(exps2)
        for iexp1 in range(0,Nexps1):
            exp1 = exps1[iexp1]
            for iexp2 in range(0,Nexps2):
                exp2 = exps2[iexp2]
                Sexps = exp2 + '-' + exp1
                print '    ' + Sexps + '...'
                owdir = cnf['ofold'] + '/' + Smods + '/' + Sexps
                sout = sub.call('mkdir -p ' + owdir, shell=True)

                # Getting the right projection in order to perform differences
                difmods = [mod1, mod2]
                difexps = [exp1, exp2]
                difreproj = [False, False]
                
                for imod in range(2):
                    mod = difmods[imod]
                    exp = difexps[imod]
                    if not gen.searchInlist(ModComProj,mod):
                        print  "  Projecting files of model '" + mod + "' exp: '" +  \
                          exp + "' with file '" + cnf['RefProjfile'] + "'"
                        reproject_modexp(mod, exp, cnf, REPRJvar, moddiffscratch,    \
                          addmoddiffs, dbg)
                        difreproj[imod] = True

                # Removing files with list of differences if should be started from 
                #   scratch or add differences
                for fdiff in difffiles:
                    difff = owdir + '/' + fdiff
                    if moddiffscratch:
                        sub.call('rm -rf' + owdir +' >& /dev/null', shell=True)
                        sub.call('rm ' + difff +' >& /dev/null', shell=True)
                        objf = open(owdir+'/all_' + fdiff,'w')
                        objf.write("## Computation and drawing of differences " +    \
                          "between '" + mod2+'@'+exp2 + "'-'" + mod1+'@'+exp1 + "'\n")
                        objf.close()
                        difff = owdir + '/all_vardiffstatistics.inf'
                        sub.call('rm ' + difff +' >& /dev/null', shell=True)
                        objf = open(owdir+'/all_' + fdiff,'w')
                        objf.write("## Computation of all differences statistics " + \
                          "between '" + mod2+'@'+exp2 + "'-'" + mod1+'@'+exp1 + "'\n")
                        objf.close()

                    if addmoddiffs:
                        if dbg:
                            print '  ' + main + ": adding model differences " +      \
                              "removing file '" + difff + "' ..."
                        sub.call('rm ' + difff +' >& /dev/null', shell=True)

                # op and var differences
                diffvarcompf = owdir + '/' + difffiles[0]
                if not os.path.isfile(owdir+'/'+difffiles[0]) or                     \
                  not os.path.isfile(owdir+'/'+difffiles[1]):
                    alldiffop, alldiffvar, Nopdiffs, Nvardiffs =                     \
                      diffvarop_listconstruct(cnf, difmods, difexps, difreproj,owdir,\
                      dbg)
                else: 
                    if dbg:
                        print warnmsg
                        print '  '+main+": getting 'op' and 'var' differences to " + \
                          "compute already from file !!"
                    alldiffop, Nopdiffops = read_diff_file(diffvarcompf)
                    alldiffvar, Nopdiffvars = read_diff_file(owdir+'/'+difffiles[1])

                # End of avoiding to repeat all the experiment search

                print "    For experiments '" + mod2+'@'+exp2 + "'-'" + mod+'@'+exp1+\
                  "' is required to compute:", Nvar, "differences"

                if dbg:
                    print '    op differences to compute _______'
                    gen.printing_dictionary(alldiffop)
                    print '    var differences to compute _______'
                    gen.printing_dictionary(alldiffvar)

# Computing differences
##
                print "    Computing operation differences ..."
                compute_op_diffs(cnf, alldiffop, owdir, moddiffscratch, dbg)
                print "    Computing variable differences ..."
                compute_var_diffs(cnf, alldiffvar, owdir, moddiffscratch, dbg)

# Plotting operation differences
##
                print "  " + main + ": Plotting operation differences' figures ..."
                dirfigf = owdir + '/diffopplotsdraw.inf'
                if figmoddiffscratch:
                    sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)

                objf = open(owdir+'/all_diffopfigures.inf','w')
                objf.write("## Drawing of all operation difference figures \n")
                objf.close()

                if addmoddifffigures:
                    if dbg:
                        print '  ' + main + ": adding model differences' figures " + \
                          "removing file '" + dirfigf + "' ..."
                    sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)

                if not os.path.isfile(dirfigf):
                    listplots, Nplt = plots_listconstruct(cnf, 'PLOTDIFFOP', dirfigf,\
                      owdir, dbg)
                else:
                    if dbg:
                        print warnmsg
                        print '  '+main+": getting plots to draw already from file !!"
                    listplots, Nplt = read_plot_file(dirfigf)

                # End of avoiding to repeat all the plots search

                print "  For experiment 'operation' differences '" + Sexps + "' is "+\
                  "required to plot:" , Nplt, "plots"

                if dbg:
                    print 'Plots to draw _______'
                    gen.printing_dictionary(listplots)

                draw_diff_plots(cnf, listplots, owdir, alldiffop, 'diffop',          \
                  figmoddiffscratch, dbg)

                # Plotting variable differences
                ##
                print "  " + main + ": Plotting variable differences' figures ..."
                dirfigf = owdir + '/diffvarplotsdraw.inf'
                if figscratch:
                    sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)

                    objf = open(owdir+'/all_diffvarfigures.inf','w')
                    objf.write("## Drawing of all variables difference figures \n")
                    objf.close()
        
                if adddifffigures:
                    if dbg:
                        print '  ' + main + ": adding differences's figures " +      \
                          "removing file '" + dirfigf + "' ..."
                    sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)

                if not os.path.isfile(dirfigf):
                    listplots, Nplt = plots_listconstruct(cnf,'PLOTDIFFVAR',dirfigf, \
                      owdir, dbg)
                else:
                    if dbg:
                        print warnmsg
                        print '  '+main+": getting plots to draw already from file !!"
                    listplots, Nplt = read_plot_file(dirfigf)

                # End of avoiding to repeat all the plots search

                print "  For experiment 'variables' differences '" + Sexps+ "' is "+ \
                  "required to plot:" , Nplt, "plots"

                if dbg:
                    print '    Plots to draw _______'
                    gen.printing_dictionary(listplots)

                draw_diff_plots(cnf, listplots, owdir, alldiffvar, 'diffvar',        \
                  figdiffscratch, dbg)

            # end of exp2 loop
        # end of exp1 loop
    # end of mod2 loop
# end of mod1 loop

### ## #
# All model/exps linear plots
### ## #

allmodexps = []
for mod in mods:
    for exp in modexps[mod]:
        allmodexps.append(mod + '/' + exp)

print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
print '  ** Multi models, multi experiments plots  '
print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
print '    Model/experiments:', allmodexps

owdir = cnf['ofold'] + '/allmodexps'
sout = sub.call('mkdir -p ' + owdir, shell=True)

# Plotting all model-experiments lines
##
print "  " + main + ": Plotting all model-experiments variable figures ..."
allf = owdir + '/allmodexp.inf'
dirfigf = owdir + '/allmodexpfigures.inf'
allfigf = owdir+'/all_figures.inf'
if figallmodexpscratch:
    sout = sub.call('rm ' + allf + ' >& /dev/null', shell=True)
    sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)
    sout = sub.call('rm ' + allfigf + ' >& /dev/null', shell=True)

    objf = open(allfigf,'w')
    objf.write("## Drawing of all variables figures for all model-experiments\n")
    objf.close()
        
if addallmodexpfigures:
    if dbg:
        print '  ' + main + ": adding model-experiment differences removing " +      \
          " file '" + allfigf + "' ..."
        print '  ' + main + ": adding model-experiment differences' figures " +      \
          "removing file '" + dirfigf + "' ..."
    sout = sub.call('rm ' + allfigf + ' >& /dev/null', shell=True)
    sout = sub.call('rm ' + dirfigf + ' >& /dev/null', shell=True)

if not os.path.isfile(dirfigf):
    listplots, Nplt = plots_listconstruct(cnf, 'PLOTALLMODEXP', dirfigf, owdir, dbg)
else:
    if dbg:
        print warnmsg
        print '  ' + main + ": getting plots to draw already from file !!"
    listplots, Nplt = read_plot_file(dirfigf)

# End of avoiding to repeat all the plots search

print "  Is required to plot:", Nplt, " all model-experiment plots"

if dbg:
    print '    Plots to draw _______'
    gen.printing_dictionary(listplots)

if not os.path.isfile(allf):
    # Constructing the information of modexp / variable
    allmodexpvar, Nallmodexpvar = allmodexps_listconstruct(cnf, allmodexps, owdir, dbg)
else:
    allmodexpvar, Nallmodexpvar = allmodexps_read(allf)

draw_allmodexp_plots(cnf, allmodexpvar, listplots, modexps, owdir, figallmodexpscratch, dbg)

print main + ': all files and figures have been properly done !!'

# Creation of a pdf with all figures
##
###
owdir = cnf['ofold']
os.chdir(owdir)

texout ='model_graphics-images'
print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
print '  ** Creation of a pdf with all figures  '
print '  ******* ****** ***** **** *** ** * ** *** **** ***** ****** *******'
print '  ' + owdir + '/' + texout + '.pdf'

if dbg:
    print main + ": generation of a pdf file '" + owdir + '/' + texout +             \
      "' with all the figures..."

otexf = open(owdir + '/' + texout + '.tex', 'w')
otexf.write('\\documentclass{article}\n')
otexf.write('\\usepackage{graphicx}\n')
otexf.write('\\usepackage[colorlinks=true,urlcolor=blue]{hyperref}\n')
otexf.write('\\textheight=23cm\n')
otexf.write('\\textwidth=18cm\n')
otexf.write('\\oddsidemargin=-1cm\n')
otexf.write('\\evensidemargin=-1cm\n')
otexf.write('\\topmargin=-1cm\n')

otexf.write('\n')
otexf.write('\\begin{document}\n')
otexf.write('\n')
otexf.write('\\def\\fontBloc{\\fontfamily{\\sfdefault}\\large\\bfseries}\n')
otexf.write('\\def\\linia{\\setlength{\\unitlength}{1mm}\n')
otexf.write('\\line(1,0){80}}\n')
otexf.write('\\newcommand{\\modexp}[1]{\\clearpage\n')
otexf.write('\\noindent\\linia\n')
otexf.write('\\section{#1}\n')
otexf.write('\\linia}\n')
otexf.write('\n')

otexf.write('\\title{model\_graphics: ' + gen.latex_text(owdir) + '}\n')
otexf.write('\\maketitle\n')
otexf.write('\n')
otexf.write('\\tableofcontents\n')
otexf.write('\\listoffigures\n')
otexf.write('\\clearpage\n')

Nmods = len(mods)
lmodexps = []
ldiffmodexps = []
# Figures from mod/exp pairs
for mod in mods:
    exps = modexps[mod]
    for exp in exps:
        print '  direct plots:', mod, exp
        lmodexps.append(mod+'/'+exp)
        owdirme = owdir + '/' + mod + '/' + exp
        otexf.write('% ' + mod + ' ' + exp + '\n')
        otexf.write('\\modexp{' + gen.latex_text(mod+' '+exp) + '}\n')
        listplots, Nplt = read_plot_file(owdirme + '/directplotsdraw.inf')
        for plot in listplots.keys():
            plots = listplots[plot]
            ops = []
            for plt in plots:
                opn = plt.split('#')[0].split('|')[1]
                if not gen.searchInlist(ops,opn): ops.append(opn)
            for op in ops:
                opS = op.replace('+','_')
                modexpimg = gen.files_folder_HMT(owdirme, head=plot, middle=opS,     \
                  tail=cnf['kindfig'])
                Nimages = len(modexpimg)
                if Nimages > 0:
                    for img in range(Nimages): modexpimg[img] = owdirme + '/' +      \
                      modexpimg[img]
                    caption = mod + ' ' + exp + ' ' + plot + ' ' + op
                    flab = 'fig:' + caption.replace(' ', '_')
                    try:
                        with gen.Capturing() as output:
                            gen.latex_fig_array(modexpimg, otexf,                    \
                              gen.latex_text(caption), flab, refsize=0.5, dist='sqr',\
                              dorest='center')
                    except:
                        print errmsg
                        print 'gen.latex_fig_array(',modexpimg,', '+otexf+', '+      \
                              'gen.latex_text(' + caption + '), ' + flab +           \
                              ", refsize=0.5, dist='sqr', dorest='center')"
                        for s1out in output: print s1out
                        quit(-1)
                    otexf.write('%\\clearpage\n')

    # differences between experiments
    Nexps = len(exps)
    for iexp1 in range(0,Nexps-1):
        exp1 = exps[iexp1]
        for iexp2 in range(iexp1+1,Nexps):
            exp2 = exps[iexp2]
            Sexps = exp2 + '-' + exp1
            owdirme = owdir + '/' + mod + '/' + Sexps
            print '    diffs:',  mod + '/' + Sexps + ' ...'
            ldiffmodexps.append(mod + '/' + Sexps)
            for diffn in ['op', 'var']:
                fileplotinf = owdirme + '/diff' + diffn + 'plotsdraw.inf'
                if os.path.isfile(fileplotinf):
                    otexf.write('% ' + mod + ' ' + Sexps + ' diff' + diffn + '\n')
                    otexf.write('\\modexp{' + gen.latex_text(mod+' '+Sexps+' diff'+  \
                      diffn)+ '}\n')
                    listplots, Nplt = read_plot_file(fileplotinf)
                    for plot in listplots.keys():
                        plots = listplots[plot]
                        ops = []
                        for plt in plots:
                            opn = plt.split('#')[0].split('|')[1]
                            if not gen.searchInlist(ops,opn): ops.append(opn)
                        for op in ops:
                            opS = op.replace('+','-')
                            modexpimg = gen.files_folder_HMT(owdirme, head=plot,     \
                              middle=opS, tail=cnf['kindfig'])
                            Nimages = len(modexpimg)
                            if Nimages > 0:
                                for img in range(Nimages): modexpimg[img] = owdirme +\
                                  '/' + modexpimg[img]
                                caption = mod + ' ' + Sexps + ' ' + plot + ' ' + op +\
                                  ' diff' + diffn 
                                flab = 'fig:' + caption.replace(' ', '_')
                                try:
                                    with gen.Capturing() as output:
                                        gen.latex_fig_array(modexpimg, otexf,        \
                                          gen.latex_text(caption), flab, refsize=0.5,\
                                          dist='sqr', dorest='center')
                                except:
                                    print errmsg
                                    print 'gen.latex_fig_array(',modexpimg,', '+     \
                                      otexf+', gen.latex_text(' + caption + '), ' +  \
                                      flab + ", refsize=0.5, dist='sqr', " +         \
                                      "dorest='center')"
                                    for s1out in output: print s1out
                                    quit(-1)

                                otexf.write('%\\clearpage\n')
print '    mods-exps diffs ...'
for imod1 in range(0,Nmods-1):
    mod1 = mods[imod1]
    exps1 = modexps[mod1]
    for imod2 in range(imod1+1,Nmods):
        mod2 = mods[imod2]
        exps2 = modexps[mod2]

        Smods = mod2 + '-' + mod1
        Nexps1 = len(exps1)
        Nexps2 = len(exps2)
        for iexp1 in range(0,Nexps1):
            exp1 = exps1[iexp1]
            for iexp2 in range(0,Nexps2):
                exp2 = exps2[iexp2]
                Sexps = exp2 + '-' + exp1
                print '      ' + Smods + ' ' + Sexps + ' ...'
                ldiffmodexps.append(Smods + '/' + Sexps)
                owdirme = owdir + '/' + Smods + '/' + Sexps
                for diffn in ['op', 'var']:
                    fileplotinf = owdirme + '/diff' + diffn + 'plotsdraw.inf'
                    if os.path.isfile(fileplotinf):
                            otexf.write('% '+Smods + ' ' + Sexps + ' diff'+diffn+'\n')
                            otexf.write('\\modexp{' + gen.latex_text(Smods+' '+Sexps+\
                              ' diff'+diffn) + '}\n')
                            listplots, Nplt = read_plot_file(fileplotinf)
                            for plot in listplots.keys():
                                plots = listplots[plot]
                                ops = []
                                for plt in plots:
                                    opn = plt.split('#')[0].split('|')[1]
                                    if not gen.searchInlist(ops,opn): ops.append(opn)
                                for op in ops:
                                    opS = op.replace('+','-')
                                    modexpimg = gen.files_folder_HMT(owdirme,        \
                                      head=plot, middle=opS, tail=cnf['kindfig'])
                                    Nimages = len(modexpimg)
                                    if Nimages > 0:
                                        for img in range(Nimages): modexpimg[img] =  \
                                          owdirme + '/' + modexpimg[img]
                                        caption = Smods + ' ' + Sexps + ' ' + plot + \
                                          ' ' + op + ' diff' + diffn 
                                        flab = 'fig:' + caption.replace(' ', '_')
                                        try:
                                            with gen.Capturing() as output:
                                                gen.latex_fig_array(modexpimg, otexf,\
                                                  gen.latex_text(caption), flab,     \
                                                  refsize=0.5, dist='sqr',           \
                                                  dorest='center')
                                        except:
                                            print errmsg
                                            print 'gen.latex_fig_array(',modexpimg,  \
                                              ', '+otexf+', gen.latex_text(' +       \
                                              caption + '), '+flab+ ", refsize=0.5,"+\
                                              " dist='sqr', dorest='center')"
                                            for s1out in output: print s1out
                                            quit(-1)
                                        otexf.write('%\\clearpage\n')

# allmodexp
owdirme = owdir + '/allmodexps'
print '    allmodexps...'
fileplotinf = owdirme + '/allmodexpfigures.inf'
if os.path.isfile(fileplotinf):
    otexf.write('% allmodexps\n')
    otexf.write('\\modexp{allmodexps}\n')
    listplots, Nplt = read_plot_file(fileplotinf)
    for plot in listplots.keys():
        plots = listplots[plot]
        ops = []
        for plt in plots:
            opn = plt.split('|')[1]
            if not gen.searchInlist(ops,opn): ops.append(opn)
        for op in ops:
            opS = op.replace('+','-')
            modexpimg = gen.files_folder_HMT(owdirme, head=plot, middle=opS,         \
              tail=cnf['kindfig'])
            Nimages = len(modexpimg)
            if Nimages > 0:
                for img in range(Nimages): modexpimg[img] = owdirme + '/' +          \
                    modexpimg[img]
                caption = 'allmodexps ' + plot + ' ' + op + ' diff' + diffn 
                flab = 'fig:' + caption.replace(' ', '_')
                try:
                    with gen.Capturing() as output:
                        gen.latex_fig_array(modexpimg,otexf, gen.latex_text(caption),\
                          flab, refsize=0.9, dist='sqr', dorest='center')
                except:
                    print errmsg
                    print 'gen.latex_fig_array(',modexpimg,', '+otexf+               \
                      ', gen.latex_text(' + caption + '), ' +flab+ ", refsize=0.9," +\
                      " dist='sqr', dorest='center')"
                    for s1out in output: print s1out
                    quit(-1)

                otexf.write('%\\clearpage\n')
#
# Grouping by type of figure
#
Nmods = len(mods)
# Figures from mod/exp pairs
owsearch = owdir + '/' + mods[0] + '/' + modexps[mods[0]][0]
print "  Looking in '" + owsearch + "' for figures by kind and variable..."
otexf.write('% figures by kind and variable from:' + owsearch + '\n')
modexpimg = gen.files_folder_HMT(owsearch, tail=cnf['kindfig'])
# dictionary with kind of plots
dplots = {}
# list with kind of plots
lplots = []
# dictionary with variable-statistics of each kind of plot
dvarstats = {} 
# list with variables of each kind of plot
lvars = []
for plot in modexpimg:
    Nvars, plotk, varn, statn = figinf_figname(plot,cnf)
    if not gen.searchInlist(lplots,plotk):
        lplots.append(plotk)
        if not dplots.has_key(plotk+'_'+varn):
            if not gen.searchInlist(lvars,varn): lvars.append(varn)
            dplots[plotk+'_'+varn] = [statn]
        else:
            vals = dplots[plotk+'_'+varn]
            dplots[plotk+'_'+varn] = vals + [statn]
    else:
        if not dplots.has_key(plotk+'_'+varn):
            if not gen.searchInlist(lvars,varn): lvars.append(varn)
            dplots[plotk+'_'+varn] = [statn]
        else:
            vals = dplots[plotk+'_'+varn]
            dplots[plotk+'_'+varn] = vals + [statn]
        
# Figures found
print '    direct figures to group:', modexpimg
print '      plots:', lplots
print '      variables:', lvars

for dplot in dplots.keys():
    plot = dplot.split('_')[0]
    var = dplot.split('_')[1]
    otexf.write('% ' + var + '\n')
    otexf.write('\\modexp{' + gen.latex_text(var) + '}\n')
    varn1 = var.split('-')[0]
    varn2 = var.split('-')[1]
    varn = varn1+'-'+varn2
    for stn12 in dplots[dplot]:
        stn1 = stn12.split('@')[0]
        stn2 = stn12.split('@')[1]
        stn = gen.lstring_values(stn12, '@',' & ')
        plotvarfigs = []
        caption = varn.replace('-',' & ') + ' ' + plot + ' ' + stn
        flab = 'fig:' + varn + '_' + plot + '_' + stn12 + '_allmodexp'
        for me in lmodexps:
            modexpn = expGC[me].label
            modn = me.split('/')[0]
            imgns = gen.files_folder_HMT(owdir+'/'+me, head=plot+'_'+varn1,          \
              middle=stn1+'-'+varn2, tail=stn2+'.'+cnf['kindfig'])
            caption = caption + ', ' + modexpn
            imgn = owdir+'/'+me+'/'+imgns[0]
            if os.path.isfile(imgn):
                plotvarfigs.append(imgn)
        if dbg:
            print '      Grouping figures:', plotvarfigs
        try:
            with gen.Capturing() as output:
                gen.latex_fig_array(plotvarfigs, otexf, gen.latex_text(caption),     \
                  flab, refsize=0.9, dist='sqr', dorest='center')
        except:
            print errmsg
            print 'gen.latex_fig_array(',plotvarfigs,', '+otexf+', gen.latex_text(' +\
               caption + '), ' + flab + ", refsize=0.9, dist='sqr', dorest='center')"
            for s1out in output: print s1out
            quit(-1)

Nmods = len(mods)
# Figures from mod/exp differences
owsearch = owdir + '/' + ldiffmodexps[0]
print "  Looking in '" + owsearch + "' for difference figures by kind and variable..."
otexf.write('% difference figures by kind and variable from:' + owsearch + '\n')
modexpimg = gen.files_folder_HMT(owsearch, tail=cnf['kindfig'])
# dictionary with kind of plots
dplots = {}
# list with kind of plots
lplots = []
# dictionary with variable-statistics of each kind of plot
dvarstats = {} 
# list with variables of each kind of plot
lvars = []
for plot in modexpimg:
    Nvars, plotk, varn, statn = figinf_figname(plot,cnf)
    if not gen.searchInlist(lplots,plotk):
        lplots.append(plotk)
        if not dplots.has_key(plotk+'_'+varn):
            if not gen.searchInlist(lvars,varn): lvars.append(varn)
            dplots[plotk+'_'+varn] = [statn]
        else:
            vals = dplots[plotk+'_'+varn]
            dplots[plotk+'_'+varn] = vals + [statn]
    else:
        if not dplots.has_key(plotk+'_'+varn):
            if not gen.searchInlist(lvars,varn): lvars.append(varn)
            dplots[plotk+'_'+varn] = [statn]
        else:
            vals = dplots[plotk+'_'+varn]
            dplots[plotk+'_'+varn] = vals + [statn]
        
# Figures found
print '    difference figures to group:', modexpimg
print '      plots:', lplots
print '      variables:', lvars

for dplot in dplots.keys():
    plot = dplot.split('_')[0]
    var = dplot.split('_')[1]
    otexf.write('% ' + var + '\n')
    otexf.write('\\modexp{diff ' + gen.latex_text(var) + '}\n')
    varn1 = var.split('-')[0]
    varn2 = var.split('-')[1]
    varn = varn1+'-'+varn2
    for stn12 in dplots[dplot]:
        stn1 = stn12.split('@')[0]
        stn2 = stn12.split('@')[1]
        stn = gen.lstring_values(stn12, '@',' & ')
        plotvarfigs = []
        caption = varn.replace('-',' & ') + ' ' + plot + ' ' + stn
        flab = 'fig:' + varn + '_' + plot + '_' + stn12 + '_alldiffmodexp'
        for me in ldiffmodexps:
            mods = me.split('/')[0]
            exps = me.split('/')[1]
            if exps.find('-') != -1:
                expn1 = exps.split('-')[0]
                expn2 = exps.split('-')[1]
            else:
                expn1 = exps
                expn2 = exps

            if mods.find('-') != -1:
                modn1 = mods.split('-')[0]
                modn2 = mods.split('-')[1]
            else:
                modn1 = mods
                modn2 = mods
            modexpn1 = expGC[modn1 + '/' + expn1].label
            modexpn2 = expGC[modn2 + '/' + expn2].label
            modexpn = modexpn1 + '-' + modexpn2

            modn = me.split('/')[0]
            imgns = gen.files_folder_HMT(owdir+'/'+me, head=plot+'_'+varn1,          \
              middle=stn1+'-'+varn2, tail=stn2+'.'+cnf['kindfig'])
            caption = caption + ', ' + modexpn
            imgn = owdir+'/'+me+'/'+imgns[0]
            if os.path.isfile(imgn):
                plotvarfigs.append(imgn)
        if dbg:
            print '      Grouping figures:', plotvarfigs
        try:
            with gen.Capturing() as output:
                gen.latex_fig_array(plotvarfigs, otexf, gen.latex_text(caption),     \
                  flab, refsize=0.9, dist='sqr', dorest='center')
        except:
            print errmsg
            print 'gen.latex_fig_array(',plotvarfigs,', '+otexf+', gen.latex_text('+ \
              caption + '), ' + flab + ", refsize=0.9, dist='sqr', dorest='center')"
            for s1out in output: print s1out
            quit(-1)

otexf.write('\\end{document}\n')
otexf.close()

try:
    with gen.Capturing() as output:
        ins = 'pdflatex ' + owdir + '/' + texout
        sout = sub.call(ins, shell=True)
except:
    print errormsg
    print '  ' + ins
    print sout
    for s1out in output: print s1out
    quit(-1)

with gen.Capturing() as output:
    sout = sub.call(ins + '>& /dev/null', shell=True)
with gen.Capturing() as output:
    sout = sub.call(ins + '>& /dev/null', shell=True)
with gen.Capturing() as output:
    sout = sub.call(ins + '>& /dev/null', shell=True)
with gen.Capturing() as output:
    sout = sub.call(ins, shell=True)
if dbg:
    print sout
    for s1out in output: print s1out

sub.call('evince ' + owdir + '/' + texout + '.pdf &', shell=True)

