# 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'

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', 'addfiles', 'addfigures', 'adddiffs',         \
      'adddifffigures', 'debug', 'models', 'ifold', 'ofold', 'warnmsg', 'errmsg',    \
      'titleoperations', 'kindfig', 'CFreftime', 'CFunitstime', 'mapval',            \
      'timekind', 'timefmt', 'timelabel' ]

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

    # Dictionary with the mandatory keys which are required as function of a given value from the dictionary
    ifkeys = {'WRFexps': ['models','WRF'], 'LMDZexps': ['models','LMDZ'],            \
      'WRF_LMDZexps': ['models','WRF_LMDZ'], 'WRFheaders': ['models','WRF'],         \
      'WRF_LMDZheaders': ['models','WRF_LMDZ'], 'LMDZheaders': ['models','LMDZ']}

    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 ExpConf[kdict]

    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 config['scratch'] == 'true':
        scr = True
        print warnmsg
        print "  " + main + ": starting from the SCRATCH !!"
        print "    10 seconds left!!"
        filescr = True
        figscr = True
        difscr = True
        figdifscr = True
        tim.sleep(10)
    else:
        scr = False
        if config['filescratch'] == 'true':
            filescr = True
            print warnmsg
            print "  " + main + ": files starting from the SCRATCH !!"
            print "    5 seconds left!!"
            tim.sleep(5)
        else:
            filescr = False

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

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

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

    if config['addfiles'] == 'true':
        addfils = True
    else:
        addfils = False

    if config['addfigures'] == 'true':
        addfigs = True
    else:
        addfigs = False

    if config['adddiffs'] == 'true':
        adddiffs = True
    else:
        adddiffs = False

    if config['adddifffigures'] == 'true':
        adddifffigs = True
    else:
        adddifffigs = False

    if config['debug'] == 'true':
        debug = True
    else:
        debug = False

    return scr, filescr, figscr, difscr, figdifscr, addfils, addfigs, adddiffs, adddifffigs, debug

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(':')
    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 variable
      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

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 sout in output: print sout
            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)

    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

    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'

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

            allvarcomp[vn + '_' + op] = [vcomp.fheader, vcomp.model, vcomp.diag, globalP]
            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() 

    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.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]]
        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 sout in output: print sout
                    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]+ "'"
                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 sout in output: print sout
                    quit(-1)
                if db:
                    for sout in output: print sout

                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':
                requiredinterpvars = ['P', 'PB', 'PSFC', 'PH', 'PHB', 'HGT', 'T',    \
                  'QVAPOR', 'XLONG', 'XLAT', 'Times']
                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 sout in output: print sout
                        quit(-1)
                    if db:
                        for sout in output: print sout

            # CFification of files
            if minf.name == 'WRF' or minf.name == 'WRF_LMDZ':
                values = vdnx + ':' + vdny + ':'+ Tref + ':' + Tunits
                try:
                    with gen.Capturing() as output:
                        ncvar.WRF_toCF(values, ifilen) 
                except:
                    print errmsg
                    print 'WRF_toCF('+ values +', ' + ifilen + ')'
                    for sout in output: print sout
                    quit(-1)
                if db:
                    for sout in output: print sout
            elif minf.name == 'LMDZ':
                ncvar.LMDZ_toCF(ifilen) 

                try:
                    with gen.Capturing() as output:
                        ncvar.LMDZ_toCF(ifilen) 
                except:
                    print errmsg
                    print 'LMDZ_toCF('+ ifilen + ')'
                    for sout in output: print sout
                    quit(-1)
                if db:
                    for sout in output: print sout
            else:
                print errmsg
                print '  ' + fname + ": no CFification for model '" +minf.name+ "' !!"
                print "    available ones: 'WRF', 'WRF_LMDZ', 'LMDZ'"
                quit(-1)

        ifile = ifile + 1

    otrackf.close()

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

        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)

    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)
        tmean: temporal mean values
        tstd: temporal standard deviation values
        turb: Taylor's turbulence decomposition value
        tvar: temporal variance values (Taylor's turbulence)
        xmean: x-axis mean 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 + "'...'"
    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 == '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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout
                
                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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout

                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':
                    values = vdnx + ':' + vdny + ':'+ Tref + ':' + Tunits
                    ncvar.WRF_toCF(values, fileon)

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

        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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout

                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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout

                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 == 'turb':
#            # turbulence values
#            vals='turbulence|'${CFvarn}
#            dims='time@time,'${dnz}'@'${vdnz}',lat@lat,lon@lon,pres@pres'
#            pyins = pyH + "/diagnostics.py -d " + dims + " -v '" +vals + \
#                "' -f " prevfiles
#            if dofile:
#                try:
#                    with gen.Capturing() as output:
#                        sout = sub.call('python ' + pyins, shell=True)
#                except:
#                    print errmsg
#                    print pyins
#                    for sout in output: print sout
#                    quit(-1)

#                if db:
#                    for sout in output: print sout

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

#            # Keeping the operations
#            pyins="python "${pyHOME}"/diagnostics.py -d "${dims}" -v '"${vals}
#            pyins=${pyins}"' -f ${cfiles}"
#            echo " " >> ${odir}/all_statsvars.inf
#            echo "# ${CFvarn}" "${vark}" >> ${odir}/all_statsvars.inf
#            echo ${pyins} >> ${odir}/all_statsvars.inf

#            varkeep=':'${CFvarn}'turb'

        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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout

                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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout

                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")

            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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout

                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")

            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 sout in output: print sout
                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 sout in output: print sout
                        quit(-1)

                    if db:
                        for sout in output: print sout

            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 sout in output: print sout
                quit(-1)

            if db:
                for sout in output: print sout

    # 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

    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]
        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)
            # Comppute variable
            compute_variable(modinf, idir, Files, odir, vinf, globalP, fscr,         \
              config['pyHOME'], config['CFreftime'], config['CFunitstime'], debug)

            # Compute variable statistics
            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

    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+'_')

    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() 

    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.split(' ')
            if values[0] == 'plots:':
                plots = gen.stringList_dictKeysVals(values[1],cV=':')
            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, 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
      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]

    if not os.path.isfile(finame):
        if db:
            print "    drawing '" + finame + "' ..."
            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 == '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]
            cline = cnt[4]

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

            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 + ')'
                print sout0
                for sout in output: print sout
                quit(-1)

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

            # keeping all figures
            trkobjf.write('\n')
            trkobjf.write("#" + 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]
            cline = cnt[4]
            # It is assumed that if the space variable is 'lon': is desired a 
            #   (lon, time) plot it 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 +';'+ cbar +  \
              ';fixsigc,' + cline + ';' + cfmt + ';' + srange + ';' + crange +',9;'+ \
              figtit + ';' + kfig + ';' + reverse + ';' + 'time|' + tunits + '|'+ \
              timekind + '|' + timefmt + '|' + timelabel

            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:
                    sout0 = sub.call(plotins, shell=True)
            except:
                print errmsg
                print drwins + '(' + graphvals + ',' + figfs + ',' + fvarS + ')'
                for sout in output: print sout
                quit(-1)

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

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

        elif '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]
            cline = cnt[4]

            # 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 = 'flip@y'

            graphvals = ','.join(CFvplot)
            graphvals = graphvals + ':' + dims + ':' + cbar + ':fixc,' + cline +     \
              ':' + cfmt + ':' + srange + ':' + crange + ',9:' + figtit + ':' +      \
              kfig + ':' + reverse + ':None'
 
            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 sout in output: print sout
                quit(-1)

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

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

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

    trkobjf.close()

    return

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']

    # 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 sout in output: print sout
                    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', '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('+','-')

                # Title of the figure:
                if ivp == 0:
                    titfigure = mod + '!' + exp + "!'" + vn + "'"
                else:
                    titfigure = titfigure + "!&!'" + vn + "'"
                opvals = op.split('+')
                for op1 in opvals:
                    if not opexplained.has_key(op1):
                        print errormsg
                        print '  '+fname+": no explanation for operation '"+op1+"' !!"
                        print '    provided:', opexplained.keys()
                        quit(-1)
                    if op1 == opvals[0]:
                        titfigure = titfigure + '$_{[' + opexplained[op1]
                        if len(opvals) == 1: titfigure = titfigure + ']}$'
                    elif op1 == opvals[len(opvals)-1]:
                        titfigure = titfigure + '\!' + opexplained[op1] + ']}$'
                    else:
                        titfigure = titfigure + '\!' + opexplained[op1]

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

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

        # End of variables-operations

    # End of kind of plots

    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' or '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

    return doopers, operationsvar, indivoperationsvar, varglobalp

def diffvarop_listconstruct(config, mods, exps, 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
      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:
        print "  'op' differences ________"
        gen.printing_dictionary(diffop)
        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]

    files1, testfiles1, allcompvar1, Nvar1= read_varcomp_file(idir1+'/varcompute.inf')
    files2, testfiles2, allcompvar2, Nvar2= read_varcomp_file(idir2+'/varcompute.inf')

    # Getting 'op' differences characteristics
    doops, opvar, indivopvar, varglobalp = get_differences_var(diffop, 'op', debug)

    diffops = {}
    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!'
                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!'
                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        

    # Getting 'var' differences characteristics
    doops, opvar, indivopvar, varglobalp = get_differences_var(diffvar, 'var', debug)

    diffvars = {}
    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!'
                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!'
                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        

    Sopdiffs = gen.dictKeysVals_stringList(diffops)
    Svardiffs = gen.dictKeysVals_stringList(diffvars)

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

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

    return diffops, diffvars, Ndiffvar, Ndiffop

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.split(' ')
            if values[0] == 'differences:':
                diffs = gen.stringList_dictKeysVals(values[1])
            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 + "'"

            dfiles = 'add|' + ifile2 + '|' + vninF + ',sub|' + ifile1 + '|' + vninF
            values = '|'.join(CFdims) + '@' + dfiles
            pyins = 'python ' + pyH + "/nc_var.py -o compute_opersvarsfiles -S '" +  \
              values + "' -v " + vninF

            try:
                with gen.Capturing() as output:
                    ncvar.compute_opersvarsfiles(values,vninF)
            except:
                print errmsg
                print 'ncvar.compute_opersvarsfiles(' + values + ',' + vninF + ')'
                for sout in output: print sout
                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()

    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
        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
        turb: Taylor's turbulence decomposition value
        tvar: temporal variance values (Taylor's turbulence)
        xmean: x-axis mean 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 == '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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout

                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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout

                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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout

                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 == 'turb':
#            # turbulence values
#            vals='turbulence|'${CFvarn}
#            dims='time@time,'${dnz}'@'${vdnz}',lat@lat,lon@lon,pres@pres'
#            pyout=`python ${pyHOME}/diagnostics.py -d ${dims} -v ${vals} -f ${cfiles}`

#            # Keeping the operations
#            pyins="python "${pyHOME}"/diagnostics.py -d "${dims}" -v '"${vals}
#            pyins=${pyins}"' -f ${cfiles}"
#            echo " " >> ${odir}/all_statsvars.inf
#            echo "# ${CFvarn}" "${vark}" >> ${odir}/all_statsvars.inf
#            echo ${pyins} >> ${odir}/all_statsvars.inf

#            varkeep=':'${CFvarn}'turb'

        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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout

                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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout

                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")

            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 sout in output: print sout
                    quit(-1)

                if db:
                    for sout in output: print sout

                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")

            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 sout in output: print sout
                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 sout in output: print sout
                        quit(-1)

                    if db:
                        for sout in output: print sout

            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 sout in output: print sout
                quit(-1)

            if db:
                for sout in output: print sout

    # End of operations
    otrackf.close()

    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 + "'"

            dfiles = 'add|' + ifile2 + '|' + vninF + ',sub|' + ifile1 + '|' + vninF
            values = '|'.join(CFdims) + '@' + dfiles
            pyins = 'python ' + pyH + "/nc_var.py -o compute_opersvarsfiles -S '" +  \
              values + "' -v " + vninF

            try:
                with gen.Capturing() as output:
                    ncvar.compute_opersvarsfiles(values,vninF)
            except:
                print errmsg
                print 'ncvar.compute_opersvarsfiles(' + values + ',' + vninF + ')'
                for sout in output: print sout
                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()

    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']

    # 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]
                mod1 = vopvals[0]
                mod2 = vopvals[1]
                mod3 = vopvals[2]
                mod4 = vopvals[3]
                headf = vopvals[4]

                modexpdiff = mod2 + '/' + exp2 + ' - ' + mod1 + '/' + exp1
                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 sout in output: print sout
                    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', '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

                if ivp == 0:
                    titfigure = modexptit + "!" + kdiff + "!'" + vn + "'"
                    vnS = vn
                else:
                    titfigure = titfigure + "!&!'" + vn + "'"
                    vnS = vnS + '!&!' + vn

                opvals = op.split('+')
                for op1 in opvals:
                    if not opexplained.has_key(op1):
                        print errormsg
                        print '  '+fname+": no explanation for operation '"+op1+"' !!"
                        print '    provided:', opexplained.keys()
                        quit(-1)
                    if op1 == opvals[0]:
                        titfigure = titfigure + '$_{[' + opexplained[op1]
                        if len(opvals) == 1: titfigure = titfigure + ']}$'
                    elif op1 == opvals[len(opvals)-1]:
                        titfigure = titfigure + '\!' + opexplained[op1] + ']}$'
                    else:
                        titfigure = titfigure + '\!' + opexplained[op1]

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

            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, odir, pyHOME,      \
              figscr, debug)

        # End of variables-operations

    # End of kind of plots

    return

# Files with information about the configuration of the script
inffiles = ['varcompute.inf', 'all_computevars.inf', 'all_statsvars.inf']

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

# 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, addfiles, addfigures, \
  adddiffs, adddifffigures, dbg = scratches(cnf)

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

# Models loop
##
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)

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

    if dbg:
        print '  model characteristics _______'
        print "  dims:", dnx, dny, dnz, dnt
        print "  var dims:", vdnx, vdny, vdnz, vdnt

    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):
            print errmsg
            print "  " + main + ": folder '" + iwdir + "' does not exist !!"
            quit(-1)

        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:
            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: 
            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:
            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:
            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:
                    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], owdir, dbg)
            else: 
                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 '"+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:
            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:
            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)
        draw_diff_plots(cnf, listplots, owdir, alldiffvar, 'diffvar', figdiffscratch,\
          dbg)

#    quit()
# end of mods loop