#!/bin/bash
# L. Fita, June 2016. Generic script to plot model outputs

main='model_graphics.bash'

# All purpouse
######## ###### ##### #### ### ## #

function list_pairs() {
# Function to provide two lists from one in order to obtain all the possibles pairs
#   listv= ':' separated list of values to construct de couples
#    $ list_pairs 1:2:3:4:5:6:7:8:9
#      list_pairs: newlists= 1|2|3|4|5|6|7|8|:2|3|4|5|6|7|8|9|
#     list_pairs: pairs= 1|2:1|3:1|4:1|5:1|6:1|7:1|8:1|9:2|3:2|4:2|5:2|6:2|7:2|8:2|9:3|4:3|5:3|6:3|7:3|8:3|9:4|5:4|6:4|7:4|8:4|9:5|6:5|7:5|8:5|9:6|7:6|8:6|9:7|8:7|9:8|9

  fname='list_pairs'

  listv=$1

  lv=`echo ${listv} | tr ':' ' '`
  Nlv=`echo ${listv} | tr ':' ' ' | wc -w | awk '{print $1}'`

  if test ${Nlv} -lt 2; then
    echo ${errormsg}
    echo "  "${fname}": too few values from list, at least 2 are required!"
    echo "    provided list: "${listv}
    exit
  fi

  Nlv_1=`expr ${Nlv} - 1`

  newlist1=`echo ${listv} | tr ':' '\n' | head -n ${Nlv_1} | tr '\n' '|'`
  newlist2=`echo ${listv} | tr ':' '\n' | tail -n ${Nlv_1} | tr '\n' '|'`

  pairs=''
  ipair=1
  il1=1
  while test ${il1} -le ${Nlv_1}; do
    lv1=`echo ${listv} | tr ':' '\n' | head -n ${il1} | tail -n 1`
    il2=`expr ${il1} + 1`
    while test ${il2} -le ${Nlv}; do
      lv2=`echo ${listv} | tr ':' '\n' | head -n ${il2} | tail -n 1`
      if test ${ipair} -eq 1; then
        pairs=${lv1}'|'${lv2}
      else
        pairs=${pairs}':'${lv1}'|'${lv2}
      fi

      ipair=`expr ${ipair} + 1`
      il2=`expr ${il2} + 1`
    done

    il1=`expr ${il1} + 1`
  done

  echo "  "${fname}": newlists= "${newlist1}':'${newlist2}
  echo "  "${fname}": pairs= "${pairs}

}

function uploadvars() {
# Function to upload variables to the system from an ASCII file as:
#   [varname] = [value]
  fileval=$1
  errormsg='ERROR -- error -- ERROR -- error'

  if test ! -f ${fileval}; then
    echo ${errormsg}
    echo "  "${fname}": file '"${fileval}"' does not exist!!"
    exit
  fi

  Nlin=`wc -l ${fileval} | awk '{print $1}'`

  ilin=1
  while test ${ilin} -le ${Nlin}; do
    line=`head -n ${ilin} ${fileval} | tail -n 1`
    varname=`echo ${line} | tr '=' ' ' | awk '{print $1}'`
    value=`echo ${line} | tr '=' ' ' | awk '{print $2}'`
    Lvarname=`expr length ${varname}'0'`

    if test ${Lvarname} -gt 1 && test ! ${varname:0:1} = '#'; then
      export ${varname}=${value}
    fi
    ilin=`expr ${ilin} + 1`
  done
}

function isInlist() {
# Function to check whether a value is in a list
  list=$1
  value=$2
  
  is=`echo ${list} | tr ':' '\n' | awk '{print "@"$1"@"}' | grep '@'${value}'@' | wc -w`
  if test ${is} -eq 1
  then
    true
  else
    false
  fi
}

function ferrmsg() {
# Function to exit and write an error message
#  comdexit: code number exit from the last execution
#  ref: script/reference from which error occurs
#  msg: message to write ('!' for spaces, '#' for end of line)
  comdexit=$1
  ref=$2
  msg=$3

  if test ${comdexit} -ne 0; then
    echo "ERROR -- error -- ERROR -- error"
    echo "  "${ref}": "$(echo ${msg} | tr '!' ' ' | tr '#' '\n')" !!"
    exit
  fi
}

function ferrmsgF() {
# Function to exit and write an error message and remove a file
#  comdexit: code number exit from the last execution
#  ref: script/reference from which error occurs
#  msg: message to write ('!' for spaces, '#' for end of line)
#  FileName: name of the file to remove in case of error
  comdexit=`expr $1 + 0`
  ref=$2
  msg=$3
  FileName=$4

  if test ${comdexit} -ne 0; then
    echo "ERROR -- error -- ERROR -- error"
    echo "  "${ref}": "$(echo ${msg} | tr '!' ' ' | tr '\#' '\n')" !!"
    rm ${FileName} >& /dev/null
    exit
  fi
}

function list_add_label(){
# Function to add the label to a list
#   list= ':' separated list to wchich add a ticket 
#   label= label to add to the list
#   ch= character to separate list values and new label
#   pos= position of the label ('beg': beginning, 'end': ending)
# $ list_add_label 1:2:3:4:5 abc # beg
# abc#1:abc#2:abc#3:abc#4:abc#5

  fname='list_add_label'

  list=$1
  label=$2
  ch=$3
  pos=$4

  ival=1
  lvs=`echo ${list} | tr ':' ' '`
  
  newlist=''
  for listv in ${lvs}; do
    if test ${pos} = 'beg'; then
      nlv=${label}${ch}${listv}
    else
      nlv=${listv}${ch}${label}
    fi
    if test ${ival} -eq 1; then
      newlist=${nlv}
    else
      newlist=${newlist}':'${nlv}
    fi
    ival=`expr ${ival} + 1`
  done

  echo ${newlist}
}

function stats_filename() {
# Function to provide the name of a given statisitcs file
#   CFvarn= CF variable name
#   fkind= kind of the file
#   headf= head of the file
# $ stats_filename ta xmean wrfout
# ta_wrfout_xmean.nc
# $ stats_filename ta pinterp@last@xmean wrfout
# ta_wrfout_pinterp_last_xmean.nc
  fname='stats_filename'

  CFvarn=$1
  fkind=$2
  headf=$3

  # Number of operations in a combo file they are '@' separated
  nopers=`echo ${fkind} | tr '@' ' ' | wc -w | awk '{print $1}'`

  if test $(index_string ${fkind} pinterp) -ne -1; then
    fhead=${headf}'p'
  else
    fhead=${headf}
  fi

  if test ${nopers} -eq 1; then
    filename=${CFvarn}'_'${fhead}'_'${fkind}
  else
    opers=`echo ${fkind} | tr '@' ' '`
    filename=${CFvarn}'_'${fhead}
    for op in ${opers}; do
      filename=${filename}'_'${op}
    done
  fi

  echo ${filename}.nc
}

function index_string(){
# Function to provide the index of a given word inside a string
#   string= string to find in
#   word= word to find
#  $ index_string 123abc456ab789bca0 bca
#  14
  fname='index_string'
  string=$1
  word=$2

  if test $# -ne 2; then
    echo ${errmsg}
    echo "  "$fname}": wrong number of arguments !!"
    echo "    2 are required and "$#" are given!"
    echo "    passed: "$0
    exit
  fi


  Lstring=`expr length ${string}`
  Lword=`expr length ${word}`
 
  if test ${Lstring} -lt 1; then
    echo ${errmsg}
    echo "  "$fname}": wrong string '"${string}"' !!"
    exit
  fi

  if test ${Lword} -lt 1; then
    echo ${errmsg}
    echo "  "$fname}": wrong word '"${word}"' !!"
    exit
  fi

  Lfinalstring=`expr ${Lstring} - ${Lword}`

  iw=0
  index=-1
  while test ${iw} -le ${Lfinalstring}; do
    if test ${string:${iw}:${Lword}} = ${word}; then
      index=${iw}
      break
    fi
    iw=`expr ${iw} + 1`
  done

  echo ${index}
  return
}

function list_filter() {
# Function to filter by a value a list
#   list: list to filter
#   value: value to filter with
#  $ list_filter '1|i:2|ii:3|iii:4|iv:5|v:1|I:2|II:3|III:4|iv:5|v' 3
#  3|iii:3|III
  fname='list_filter'

  list=$1
  varn=$2

  lvs=`echo ${list} | tr ':' ' '`

  newlist=''
  il=1
  for lv in ${lvs}; do
    if test $(index_string ${lv} ${varn}) -ne -1; then
      if test ${il} -eq 1; then
        newlist=${lv}
      else
        newlist=${newlist}':'${lv}
      fi
      il=`expr ${il} + 1`
    fi
  done

  echo ${newlist}

  return
}

function join_lists() {
# Function to join lists without repeating values
#   list1: first list to join
#   list2: second list to join
#  $ join_lists 1:2:3:4:5:6:7:8 a:b:3:c:5:6:7:9
#  1:2:3:4:5:6:7:8:a:b:c:9

  fname='join_lists'
  list1=$1
  list2=$2

  ls1s=`echo ${list1} | tr ':' ' '`
  ls2s=`echo ${list2} | tr ':' ' '`

  newlist=${list1}
  for ls2 in ${ls2s}; do
    if ! $(isInlist ${list1} ${ls2}); then
      newlist=${newlist}':'${ls2}
    fi
  done

  echo ${newlist}
}

# Specific
######## ###### ##### #### ### ## #

function variable_compute(){
# 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= '@' separated list of '[head]|[file] test files
  fname='variable_compute'

  iwdir=$1
  var=$2
  filtests=$3

  ftsts=`echo ${filtests} |tr '@' ' '`
  cancompute=true

  for ftest in ${ftsts}; do
    headerf=`echo ${ftest} | tr '|' ' ' | awk '{print $1}'`
    filen=`echo ${ftest} | tr '|' ' ' | awk '{print $2}'`
    compute=`${pyHOME}/nc_var.py -o computevar_model -f ${iwdir}/${filen} -S ${var}`
    ferrmsg $? ${main} "python!'computevar_model'!failed!#"$(echo ${compute} |       \
      tr ' ' '!')
    varmod=`echo ${compute} | tr ' ' '#' | tr '|' ' ' | awk '{print $2}' |           \
      tr '@' ' ' | awk '{print $2}' | tr '=' ' ' | awk '{print $2}'`
    vardiag=`echo ${compute} | tr ' ' '#' | tr '|' ' ' | awk '{print $2}' |          \
      tr '@' ' ' | awk '{print $3}' | tr '=' ' ' | awk '{print $2}'`
    echo "  "${var}" mod:"${varmod}" diag: "${vardiag}
    if test ${varmod} = 'None' && test ${vardiag} = 'None'; then 
      cancompute=false
    else
      cancompute=true
    # Should be considered that variable can also be computed by both ways?
      break
    fi
  done
  if ! ${cancompute}; then
    msg="there!is!no!way!to!compute!'"${var}"'!for!model!"${mod}
# Too extrict!
#    ferrmsg 1 ${main} ${msg}
    echo ${warnmsg}
    echo $(echo $msg | tr '!' ' ')
  fi

  # 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
  varcomp=${var}'|'${vark}'|'${headerf}'|'${varmod}'|'${vardiag}

  echo "varcomp= "${varcomp}
}

function compvars_listconstruct() {
# Function to construct the list of variables to compute
#   VarKinds= ':' separated list of kind variables to compute
#   TestFiles= ':' list of files to test how to compute the variable
##
  fname='compvars_listconstruct'

  VarKinds=$1
  TestFiles=$2

  varks=`echo ${VarKinds} | tr ':' ' '`

  combo=false
  ik=1
  itotv=1
  for vark in ${varks}; do
    if test ! ${vark} = 'diff'; then
      echo "    "${vark}" ..."
      case ${vark} in
        'acc')
          varvks=${varacc}
        ;;
        'combo')
          combo=true
          varvcombo=${varcombo}
          varcos=`echo ${varvcombo} | tr ':' ' '`
          varvks=''
          ivco=1
          for vco in ${varcos}; do
            vn=`echo ${vco} | tr '@' ' ' | awk '{print $1}'`
            if test ${ivco} -eq 1; then
              varvks=${vco}
            else
              varvks=${varvks}':'${vco}
            fi
            ivco=`expr ${ivco} + 1`
          done
        ;;
        'direct')
          varvks=${vardirect}
        ;;
        'last')
          varvks=${varlast}
        ;;
        'Lmean')
          varvks=${varLmean}
        ;;
        'Lsec')
          varvks=${varLsec}
        ;;
        'lmean')
          varvks=${varlmean}
        ;;
        'lsec')
          varvks=${varlsec}
        ;;
        'pinterp')
          varvks=${varpinterp}
        ;;
        'tmean')
          varvks=${vartmean}
        ;;
        'xmean')
          varvks=${varxmean}
        ;;
        'ymean')
          varvks=${varymean}
        ;;
        'zsum')
          varvks=${varzsum}
        ;;
        '*')
          echo ${errmsg}
          echo "  "${main}": variable kind '"${vark}"' not ready!!"
          exit
        ;;
      esac

# Do we have variables for this kind?
      Lvarvks=`expr length ${varvks}'0'`
      if test ${Lvarvks} -lt 2; then
        ferrmsg 1 ${main} "variable!kind!"${vark}"!without!variables"
      fi
      if test ${ik} -eq 1; then
        allvars=${varvks}
      else
        allvars=${allvars}':'${varvks}
      fi
      ik=`expr ${ik} + 1`
      vars=`echo ${varvks} | tr ':' ' '`

# Variables loop
##
      iv=1
      for var in ${vars}; do
        echo "      "${var}
        # How to compute it?
        varcomp=`variable_compute ${iwdir} ${var} ${TestFiles} | grep varcomp |  \
          awk '{print $2}'`
  
        if test ${itotv} -eq 1; then
          varcompute=${varcomp}
        else
          varcompute=${varcompute}';'${varcomp}
        fi

        iv=`expr ${iv} + 1`
        itotv=`expr ${itotv} + 1`

# end loop of variables
      done
    fi
# end of kind of variables
  done

  echo ${fname}" varcompute= "${varcompute}_${itotv}_${combo}
}

function WRF_toCF() {
# Function to pass a WRF original file to CF-conventions
#   wrff= WRF file
#   wrfvl= WRF longitude var name (it might come from WPS!)
#   wrfvL= WRF latitude var name (it might come from WPS!)

  wrff=$1
  wrfvl=$2
  wrfvL=$3

  fdims=`${pyHOME}/nc_var.py -o idims -f ${wrff} | grep alldims | tr '=' ' ' |       \
    awk '{print $3}'`

  CFdims='west_east@lon:south_north@lat:Time@time'
  CFds=`echo ${CFdims} | tr ':' ' '`
  
  pyout=`${pyHOME}/nc_var.py -o WRF_CFtime_creation -S 19491201000000,minutes        \
    -f ${wrff} -v time`
  pyn=$?
  Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
  ferrmsgF ${pyn} ${fname} "python!'WRF_CFtime_creation'!failed"${Spyout} ${wrff}
  
  pyout=`${pyHOME}/nc_var.py -o WRF_CFlonlat_creation -v ${wrfvl},${wrfvL}           \
    -f ${wrff} -S lon,lat`
  pyn=$?
  Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
  ferrmsgF ${pyn} ${fname} "python!'WRF_CFlonlat_creation'!failed"${Spyout} ${wrff}

  for CFd in ${CFds}; do
    cd=`echo ${CFd} | tr '@' ' ' | awk '{print $2}'`
    wd=`echo ${CFd} | tr '@' ' ' | awk '{print $1}'`
    if $(isInlist ${fdims} ${wd}); then
      pyout=`${pyHOME}/nc_var.py -o chdimname -f ${wrff} -S ${wd}':'${cd}`
      pyn=$?
      Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
      ferrmsgF ${pyn} ${fname} "python!'chdimname'!'"${wd}"'!failed#"${Spyout} ${wrff}
    fi
  done

  # CF attributes
  dimvars='lon:lat:time:Times'

  allfilevars=`python ${pyHOME}/nc_var.py -o ivars -f ${wrff} | grep allvars | awk '{print $3}'`
  allvs=`echo ${allfilevars} | tr ':' ' '`
  for vn in ${allvs}; do
    if ! $(isInlist ${dimvars} ${vn}); then
      varattrs=`python ${pyHOME}/generic.py -o variables_values -S ${vn}`
      stn=`echo ${varattrs} | tr ':' ' ' | awk '{print $2}'`
      lon=`echo ${varattrs} | tr ':' ' ' | awk '{print $5}' | tr '|' '!'`
      un=`echo ${varattrs} | tr ':' ' ' | awk '{print $6}' | tr '|' '!'`

      pyout=`python ${pyHOME}/nc_var.py -o varaddattr -f ${wrff} -S 'standard_name|'${stn} -v ${vn}`
      pyn=$?
      Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
      ferrmsgF ${pyn} ${fname} "python!'varaddattr'!'standard_name'!failed#"${Spyout} ${wrff}
      pyout=`python ${pyHOME}/nc_var.py -o varaddattr -f ${wrff} -S 'long_name|'${lon} -v ${vn}`
      pyn=$?
      Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
      ferrmsgF ${pyn} ${fname} "python!'varaddattr'!'long_name'!failed#"${Spyout} ${wrff}
      pyout=`python ${pyHOME}/nc_var.py -o varaddattr -f ${wrff} -S 'units|'${un} -v ${vn}`
      pyn=$?
      Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
      ferrmsgF ${pyn} ${fname} "python!'varaddattr'!'units'!failed#"${Spyout} ${wrff}
    fi
  done
}

function cleaning_varsfile() {
# Function to keep only a list of variables from a file
#   filen= file to clean
#   keepvars= ':' separated list of variables to keep
  fname='cleaning_varsfile'

  filen=$1
  keepvars=$2

  fvars=`${pyHOME}/nc_var.py -o ivars -f ${filen} | grep allvars | tr '=' ' ' |      \
    awk '{print $3}'`
  fvs=`echo ${fvars} | tr ':' ' '`

  for fv in ${fvs}; do
    if ! $(isInlist ${keepvars} ${fv}); then
      pyout=`python ${pyHOME}/nc_var.py -o varrm -v ${fv} -f ${filen}`
      pyn=$?
      Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
      ferrmsg ${pyn} ${fname} "python!'varrm'!'${fv}'!failed#"${Spyout}
    fi
  done
}

function compute_variable() {
# Function to compute a variable
#   idir: directory with the input files
#   usefiles: ',' list of files as [headerf]@[file1]|...|[fileN]
#   odir: directory to write the output files
#   cvar: [CFvarn]|[varkind]|[headerf]|[varmod]|[vardiag]
#     [CFvarn]: CF name of the variable
#     [vark]: kind of variable: 
#        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]) 
#        lsec: longitudinal section (longitudinal value must be given, [var]@[lat])
#        lmean: longitudinal mean values
#        pinterp: pressure interpolation (to the given $plevels)
#        tmean: temporal mean values
#        tmean: temporal mean values
#        xmean: x-axis mean values
#        ymean: y-axis mean values
#        zsum: vertical aggregated values
#     [headerf]: header of the files to use
#     [modvar]: variable to use from the file
#     [diagvar]: variable computed from the diagnostics (diagnostic.py)
#   mdim: name of the dimensions in the model
#   mvdim: name of the vaariable-dimensions in the model
#   scratch: whether is has to be done from the scratch or not
  fname='compute_variable'
  idir=$1
  usefiles=$2
  odir=$3
  cvar=$4
  mdim=$5
  mvdim=$6
  scratch=$7

  CFvarn=`echo ${cvar} | tr '|' ' ' | awk '{print $1}'`
  vark=`echo ${cvar} | tr '|' ' ' | awk '{print $2}'`
  headerf=`echo ${cvar} | tr '|' ' ' | awk '{print $3}'`
  modvar=`echo ${cvar} | tr '|' ' ' | awk '{print $4}'`
  diagvar=`echo ${cvar} | tr '|' ' ' | awk '{print $5}'`

  cfiles=`echo ${usefiles} | tr ',' '\n' | grep ${headerf} | tr '|' ' ' | \
    awk '{print $2}' | tr '@' ' '`

# dimensions
  dnx=`echo ${mdim} | tr ',' ' ' | awk '{print $1}'`
  dny=`echo ${mdim} | tr ',' ' ' | awk '{print $2}'`
# var-dimensions
  vdnx=`echo ${mvdim} | tr ',' ' ' | awk '{print $1}'`
  vdny=`echo ${mvdim} | tr ',' ' ' | awk '{print $2}'`

  cd ${odir}

  # Computing separately and then joinging for all files
  Ntotfiles=`echo ${cfiles} | wc -w | awk '{print $1}'`
  ifile=1

  for cf in ${cfiles}; do
    ifS=`printf "%0${Ntotfiles}d" ${ifile}`

# Computing variable
    # Changing file head when it is a pressure-interpolated variable
    if test ${vark} == 'pinterp'; then
      fhead=${headerf}p
    else
      fhead=${headerf}
    fi

    filen=${odir}/${CFvarn}_${fhead}_${ifile}-${Ntotfiles}.nc
    if ${scratch}; then 
      rm ${filen} >& /dev/null
      rm ${odir}/${CFvarn}_${fhead}.nc >& /dev/null
    fi

    if test ! -f ${filen} && test ! -f ${odir}/${CFvarn}_${fhead}.nc; then
# 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 test ! ${modvar} = 'None' && test ${diagvar} = 'None'; then
        # model variable
        values=${modvar}',0,-1,-1'
        vs=${modvar},${vdnx},${vdny},${vdnz},${vdnt}
        pyout=`${pyHOME}/nc_var.py -f ${cf} -o DataSetSection_multivars -v ${vs}     \
          -S ${values} | grep succesfull | awk '{print $6}' | tr '"' ' '`
        pyn=$?
        Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
        ferrmsg ${pyn} ${fname} "python!'DataSetSection_multivars'!failed"${Spyout}
        mv ${pyout} ${filen}

        # Keeping the operations
        pyins=${pyHOME}"/nc_var.py -f "${cf}" -o DataSetSection_multivars -v "${vs}
        pyins=${pyins}" -S "${values}
        echo " " >> ${odir}/all_computevars.inf
        echo "# ${CFvarn}" "${modvar}" >> ${odir}/all_computevars.inf
        echo ${pyins} >> ${odir}/all_computevars.inf

        pyout=`${pyHOME}/nc_var.py -f ${filen} -o chvarname -v ${modvar} -S ${CFvarn}`
        pyn=$?
        Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
        ferrmsgF ${pyn} ${fname} "python!'chvarname'!failed"${Spyout} ${filen}
      else
        # diagnostic variable
        dims=${dnt}@${vdnt},${dnz}@${vdnz},${dny}@${vdny},${dnx}@${vdnx}
        diagn=`echo ${diagvar} | tr ':' '\n' | head -n 1`
        Ndiagvars=`echo ${diagvar} | tr ':' ' ' | wc -w | awk '{print $1}'`
        idiagv=2
        diagc=''
        while test ${idiagv} -le ${Ndiagvars}; do
          diag1=`echo ${diagvar} | tr ':' '\n' | head -n ${idiagv} | tail -n 1`
          if test ${idiagv} -eq 2; then
            diagc=${diag1}
          else
            diagc=${diagc}'@'${diag1}
          fi
          idiagv=`expr ${idiagv} + 1`
        done
        pyout=`python ${pyHOME}/diagnostics.py -f ${cf} -d ${dims} -v ${diagn}'|'${diagc}`
        pyn=$?
        Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
        ferrmsg ${pyn} ${fname} "python!'diagnostics'!failed#"${Spyout}
        mv diagnostics.nc ${filen}

        # Keeping the operations
        pyins=${pyHOME}"/diagnostics.py -f "${cf}" -d "${dims}" -v '"${diagn}"|"
        pyins=${pyins}${diagc}"'"
        echo " " >> ${odir}/all_computevars.inf
        echo "# ${CFvarn}" "${diagvar}" >> ${odir}/all_computevars.inf
        echo ${pyins} >> ${odir}/all_computevars.inf
      fi

      # adding CF lon,lat,time in WRF files
      if test ${headerf:0:3} = 'wrf'; then
        WRF_toCF ${filen} ${vdnx} ${vdny}
      fi
      # Attaching necessary variables for the pressure interpolation
      if test ${vark} == 'pinterp'; then
        requiredinterpvars='P:PB:PSFC:PH:PHB:HGT:T:QVAPOR:'
        rqvs=`echo ${requiredinterpvars} | tr ':' ' '`
        echo "  "${fname}": adding variables: "${rqvs}" to allow pressure interpolation"
        for rqv in ${rqvs}; do
          pyout=`${pyHOME}/nc_var.py -o fvaradd -S ${cf},${rqv} -f ${filen}`
          pyn=$?
          Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
          ferrmsgF ${pyn} ${fname} "python!'fvaradd'!failed#"${Spyout} ${filen}
        done
      fi
    fi

    ifile=`expr ${ifile} + 1`
# End of files
  done

  # Joining variable files
  filen=${odir}/${CFvarn}_${fhead}.nc
  if ${scratch}; then rm ${filen} >& /dev/null; fi

  if test ! -f ${filen}; then
    pyout=`python ${pyHOME}/nc_var.py -f ${CFvarn}'_'${fhead}'_,-,.nc'               \
      -o netcdf_fold_concatenation_HMT -S ./,time -v all`
    pyn=$?
    Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
    ferrmsg ${pyn} ${fname} "python!'netcdf_fold_concatenation_HMT'!failed#"${Spyout}
    mv netcdf_fold_concatenated_HMT.nc ${filen}
    if test -f ${filen}; then
      rm ${CFvarn}_${fhead}_*-*.nc
    fi
  fi
}

function compute_vars() {
# Function to compute variables
#   Iwdir= input folder
#   Owdir= output folder
#   Files= files to use
#   Moddims= list of names of dimensions of the file
#   Modvdims= list of names of the dimensions variables of the file
#   VarCompute= list of variables to compute as ';' list of 
#     [CFvarn]|[varkind]|[headerf]|[varmod]|[vardiag]
#   Filescratch= wether files have to be compute from the scratch or not
    fname='compute_vars'

    Iwdir=$1
    Owdir=$2
    Files=$3
    Moddims=$4
    Modvdims=$5
    VarCompute=$6
    Filescratch=$7

    cd $Owdir
    ic=1
    isc=1
    cvars=`echo ${VarCompute} | tr ';' ' '`
    for cvar in ${cvars}; do
      CFv=`echo ${cvar} | tr '|' ' ' | awk '{print $1}'`
      vark=`echo ${cvar} | tr '|' ' ' | awk '{print $2}'`
      fileh=`echo ${cvar} | tr '|' ' ' | awk '{print $3}'`
      modv=`echo ${cvar} | tr '|' ' ' | awk '{print $4}'`
      diagv=`echo ${cvar} | tr '|' ' ' | awk '{print $5}'`
      if ${dbg}; then echo "        "${CFv}"; "${modv}" "${diagv}; fi

      if test ! ${modv} = 'None' || test ! ${diagv} = 'None'; then
        compute_variable ${Iwdir} ${Files} ${Owdir} ${cvar} ${Moddims} ${Modvdims}   \
          ${Filescratch}
        ic=`expr ${ic} + 1`

        if test ! ${vark} = 'diff' && test ! ${vark} = 'combo'; then
          if test ${vark} = 'pinterp'; then
            fhead=${fileh}'p'
          else
            fhead=${fileh}
          fi
          compute_statistics ${iwdir} ${CFv}_${fhead}.nc ${owdir} ${cvar}            \
            ${moddims} ${modvdims} ${filescratch}
          isc=`expr ${isc} + 1`
        else
          if test ${vark} = 'diff'; then
            echo "  "${main}": differences will be calculated when all the "         \
              "model/experiments will be done !"
          elif test ${vark} = 'combo'; then
            echo "  "${main}": combos will be calculated later when all the "        \
              "variables will be done !"
          fi
        fi
      else
        if ${dbg}; then echo "        not '"${CFv}"' for model '"${mod}"' !!"; fi
      fi

      # exit
# end of computing vars
    done
}

function compute_statistics(){
# Function to compute different statistics
#   idir: directory with the input files
#   usefiles: ',' list of files to use [file1],...,[fileN]
#   odir: directory to write the output files
#   cvar: [CFvarn]|[varkind]|[headerf]|[varmod]|[vardiag]
#     [CFvarn]: CF name of the variable
#     [vark]: kind of variable: 
#        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
#        turb: Taylor's turbulence decomposition value
#        xmean: x-axis mean values
#        ymean: y-axis mean values
#        zsum: vertical aggregated values
#     [headerf]: header of the files to use
#     [modvar]: variable to use from the file
#     [diagvar]: variable computed from the diagnostics (diagnostic.py)
#   mdim: name of the dimensions in the model
#   mvdim: name of the vaariable-dimensions in the model
#   scratch: whether is has to be done from the scratch or not
  fname='compute_statistics'

  idir=$1
  usefiles=$2
  odir=$3
  cvar=$4
  mdim=$5
  mvdim=$6
  scratch=$7

  # Getting a previous file name to continue with(for combinations)
  if test $# -eq 8; then
    echo ${warnmsg}
    echo "  "${fname}": using a previous file '"$8"' to continue with!!"
    prevfile=$8
    usefiles=$8
  fi

  CFvarn=`echo ${cvar} | tr '|' ' ' | awk '{print $1}'`
  vark=`echo ${cvar} | tr '|' ' ' | awk '{print $2}'`
  headerf=`echo ${cvar} | tr '|' ' ' | awk '{print $3}'`
  modvar=`echo ${cvar} | tr '|' ' ' | awk '{print $4}'`
  diagvar=`echo ${cvar} | tr '|' ' ' | awk '{print $5}'`

  cfiles=`echo ${usefiles} | tr ',' ' '`

# dimensions
  dnx=`echo ${mdim} | tr ',' ' ' | awk '{print $1}'`
  dny=`echo ${mdim} | tr ',' ' ' | awk '{print $2}'`
# var-dimensions
  vdnx=`echo ${mvdim} | tr ',' ' ' | awk '{print $1}'`
  vdny=`echo ${mvdim} | tr ',' ' ' | awk '{print $2}'`

  cd ${odir}

  if test $# -ne 8; then
    if test ${vark} == 'pinterp'; then
      fhead=${headerf}p
    else
      fhead=${headerf}
    fi
    filen=${CFvarn}_${fhead}_${vark}.nc
  else
    Lprevfile=`expr length ${prevfile}`
    Lprevfile3=`expr ${Lprevfile} - 3`
    filen=${prevfile:0:$Lprevfile3}_${vark}.nc
  fi
  if ${scratch}; then rm ${filen} >& /dev/null; fi

  if test ! -f ${filen}; then
    # Computing stats
    case ${vark} in
       # temporal accumulated values
      'acc') 
        echo "  "${fname}": kind '"${vark}"' not ready !!"
        exit
      ;;
      # differences between models
      'diff')
        echo "  "${fname}": kind '"${vark}"' not ready !!"
        exit
      ;;
      # no statistics
      'direct')
        cp ${cfiles} ${filen}
      ;;
      # last temporal value
      'last')
        Lcfiles=`expr length ${cfiles}`
        Lcfiles3=`expr ${Lcfiles} - 3`
        vals='time,-9,0,0'
        pyout=`python ${pyHOME}/nc_var.py -o DataSetSection -S ${vals} -f ${cfiles}  \
          -v ${CFvarn}`
        pyn=$?
        Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
        ferrmsg ${pyn} ${fname} "python!'DataSetSection'!'last'!failed#"${Spyout}
        mv ${cfiles:0:${Lcfiles3}}_time_B-9-E0-I0.nc ${filen}

        # Keeping the operations
        pyins=${pyHOME}"/nc_var.py -o DataSetSection -S "${vals}" -f "${cfiles}
        pyins=${pyins}"-v "${CFvarn}
        echo " " >> ${odir}/all_statsvars.inf
        echo "# ${CFvarn}" "${vark}" >> ${odir}/all_statsvars.inf
        echo ${pyins} >> ${odir}/all_statsvars.inf
      ;;
      # latitudinal mean values
      'Lmean')
        echo "  "${fname}": kind '"${vark}"' not ready !!"
        exit
      ;;
      # latitudinal section (latitudinal value must be given, [var]@[lat]) 
      'Lsec')
        echo "  "${fname}": kind '"${vark}"' not ready !!"
        exit
      ;;
      # longitudinal section (longitudinal value must be given, [var]@[lat])
      'lsec')
        echo "  "${fname}": kind '"${vark}"' not ready !!"
        exit
      ;;
      # longitudinal mean values
      'lmean')
        echo "  "${fname}": kind '"${vark}"' not ready !!"
        exit
      ;;
      # pinterp: pressure interpolation (to the given $plevels)
      'pinterp')
        vals=${plevels}',1,1'
        pyout=`python $pyHOME/nc_var.py -o pinterp -f ${cfiles} -S ${vals}           \
          -v ${CFvarn}`
        pyn=$?
        Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
        ferrmsg ${pyn} ${fname} "python!'pinterp'!failed#"${Spyout}
        cp pinterp.nc ${filen}

        # Keeping the operations
        pyins=${pyHOME}"/nc_var.py -o pinterp -S "${vals}" -f "${cfiles}
        pyins=${pyins}"-v "${CFvarn}
        echo " " >> ${odir}/all_statsvars.inf
        echo "# ${CFvarn}" "${vark}" >> ${odir}/all_statsvars.inf
        echo ${pyins} >> ${odir}/all_statsvars.inf

        # adding CF lon,lat,time in WRF files
        if test ${headerf:0:3} = 'wrf'; then
          WRF_toCF ${filen} ${vdnx} ${vdny}
        fi
      ;;
      # temporal mean values
      'tmean')
        vals='time|-1,time,mean,lon:lat:'${vdnz}':time:pres'
        dims='time@time,'${dnz}'@'${vdnz}',lat@lat,lon@lon'

        pyout=`python ${pyHOME}/nc_var.py -o file_oper_alongdims -S ${vals}          \
          -f ${cfiles} -v ${CFvarn}`
        pyn=$?
        Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
        ferrmsg ${pyn} ${fname} "python!'file_oper_alongdims'!'tmean'!failed#"${Spyout}
        mv file_oper_alongdims_mean.nc ${filen}

        # Keeping the operations
        pyins=${pyHOME}"/nc_var.py -o file_oper_alongdims -S '"${vals}"' -f "${cfiles}
        pyins=${pyins}" -v "${CFvarn}
        echo " " >> ${odir}/all_statsvars.inf
        echo "# ${CFvarn}" "${vark}" >> ${odir}/all_statsvars.inf
        echo ${pyins} >> ${odir}/all_statsvars.inf

        varkeep=':'${CFvarn}'mean:timestats'
      ;;
      # turbulence values
      'turb')
        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}`
        pyn=$?
        Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
        ferrmsg ${pyn} ${fname} "python!'diagnostics'!'turbulence'!failed#"${Spyout}
        mv diagnostics.nc ${filen}

        # 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'
      ;;
      # x-axis mean values
      'xmean')
        vals='lon|-1,lon,mean,lon:lat:'${vdnz}':time:pres'
        dims='time@time,'${dnz}'@'${vdnz}',lat@lat,lon@lon'

        pyout=`python ${pyHOME}/nc_var.py -o file_oper_alongdims -S ${vals}         \
          -f ${cfiles} -v ${CFvarn}`
        pyn=$?
        Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
        ferrmsg ${pyn} ${fname} "python!'file_oper_alongdims'!'tmean'!failed#"${Spyout}
        mv file_oper_alongdims_mean.nc ${filen}

        # Keeping the operations
        pyins=${pyHOME}"/nc_var.py -o file_oper_alongdims -S '"${vals}"' -f "${cfiles}
        pyins=${pyins}" -v "${CFvarn}
        echo " " >> ${odir}/all_statsvars.inf
        echo "# ${CFvarn}" "${vark}" >> ${odir}/all_statsvars.inf
        echo ${pyins} >> ${odir}/all_statsvars.inf

        varkeep=':'${CFvarn}'mean:lonstats'
      ;;
      # y-axis mean values
      'ymean')
        vals='lat|-1,lat,mean,lon:lat:'${vdnz}':time:pres'
        dims='time@time,'${dnz}'@'${vdnz}',lat@lat,lon@lon'

        pyout=`python ${pyHOME}/nc_var.py -o file_oper_alongdims -S ${vals}         \
          -f ${cfiles} -v ${CFvarn}`
        pyn=$?
        Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
        ferrmsg ${pyn} ${fname} "python!'file_oper_alongdims'!'tmean'!failed#"${Spyout}
        mv file_oper_alongdims_mean.nc ${filen}

        # Keeping the operations
        pyins=${pyHOME}"/nc_var.py -o file_oper_alongdims -S '"${vals}"' -f "${cfiles}
        pyins=${pyins}" -v "${CFvarn}
        echo " " >> ${odir}/all_statsvars.inf
        echo "# ${CFvarn}" "${vark}" >> ${odir}/all_statsvars.inf
        echo ${pyins} >> ${odir}/all_statsvars.inf

        varkeep=':'${CFvarn}'mean:latstats'
      ;;
      # vertical aggregated values
      'zsum')
        echo "  "${fname}": kind '"${vark}"' not ready !!"
        exit
      ;;
      *)
        echo ${errmsg}
        echo "  "${fname}": kind '"${vark}"' not ready !!"
        exit
      ;;
    esac
    cleaning_varsfile ${filen} ${CFvarn}':lon:lat:pres:time'${varkeep}
  fi

}

function compute_combos(){
# Function to compute the given list of combos
#   Iwdir= input folder
#   Owdir= output folder
#   Files= files to use
#   Testfiles= files for testing the wat of compute a given variable
#   Moddims= dimensions of the input files
#   Modvdims= variable dimensions of the input files
#   VarCombo= ':' separated list of combos to compute
#   Filescratch= do we need to start files from the scratch?
  fname='compute_combos'

  Iwdir=$1
  Owdir=$2
  Files=$3
  Testfiles=$4
  Moddims=$5
  Modvdims=$6
  VarCombo=$7
  Filescratch=$8

  varcos=`echo ${VarCombo} | tr ':' ' '`
  for vco in ${varcos}; do
    echo "        "${vco}
    CFvarn=`echo ${vco} | tr ';' ' ' | awk '{print $1}'`
    cvar=`variable_compute ${Iwdir} ${CFvarn} ${Testfiles} | grep varcomp |          \
      awk '{print $2}'`
    CFv=`echo ${cvar} | tr '|' ' ' | awk '{print $1}'`
    vark=`echo ${cvar} | tr '|' ' ' | awk '{print $2}'`
    fileh=`echo ${cvar} | tr '|' ' ' | awk '{print $3}'`
    modv=`echo ${cvar} | tr '|' ' ' | awk '{print $4}'`
    diagv=`echo ${cvar} | tr '|' ' ' | awk '{print $5}'`

#        if $(index_string ${vark} 'pinterp') -ne -1; then
#          fhead=${fileh}'p'
#        else
#          fhead=${fileh}
#        fi
    combs=`echo ${vco} | tr ';' ' ' | awk '{print $2}' | tr '@' ' '`
    ifile=${CFv}_${fhead}
    for comb in ${combs}; do
      combcvar=${CFv}'|'${comb}'|'${fhead}'|'${modv}'|'${diagv}
      compute_statistics ${Iwdir} ${ifile}.nc ${Owdir} ${combcvar} ${Moddims}    \
        ${Modvdims} ${Filescratch} ${ifile}.nc
      ifile=${ifile}_${comb}

      # Adding 'surnames' to the variable's name
      if test ${comb} = 'turb'; then
        CFv=${CFv}'turb'
      fi
      if $(isInlist ${varmeanname} ${comb}); then
        CFv=${CFv}'mean'
      fi

      isc=`expr ${isc} + 1`
    # End of combo 
    done
  # End of combo variables
  done

  return
}

function createlist_plots(){
# Function to create the list of plots to draw
#   Drawplots= ':' separated list of plots to draw
#   KindPlots= ':' list of plots to draw this time

  fname='createlist_plots'

  Drawplots=$1
  KindPlots=$2

  drwplts=`echo ${Drawplots} | tr ':' ' '`
  ikp=1
  plots=''
  for drw in ${drwplts}; do
    if $(isInlist ${KindPlots} ${drw}); then
      case ${drw} in
        'map2Dsfc')
          plts=${pltmap2Dsfc}
        ;;
        'map3D')
          plts=${pltmap3D}
        ;;
        'shadcont2Dsfc')
          plts=${pltshadcont2Dsfc}
        ;;
        'shadconthovmsfc')
          plts=${pltshadconthovmsfc}
        ;;
        'shadcont2Dzsec')
          plts=${pltshadcont2Dzsec}
        ;;
        *)
          echo ${errmsg}
          echo "  "${main}": plot kind '"${drw}"' not ready !!"
          exit
      esac
      # Do we have plots for this kind?
      Lkplots=`expr length ${plts}'0'`
      if test ${Lkplots} -lt 2; then
        ferrmsg 1 ${main} "plot!kind!"${drw}"!without!variables"
      fi
      kplots=`list_add_label ${plts} ${drw} + beg`
      if test ${ikp} -eq 1; then
        plots=${kplots}
      else
        plots=${plots}':'${kplots}
      fi
      ikp=`expr ${ikp} + 1`
    fi
# End of coincident direct plots
  done

  echo ${fname}" listplots= "${plots}_${ikp}
}

function createlist_diffplots(){
# Function to create the list of diff plots to draw
#   Drawplots= ':' separated list of plots to draw
#   KindPlots= ':' list of plots to draw this time

  fname='createlist_diffplots'

  Drawplots=$1
  KindPlots=$2

  drwplts=`echo ${Drawplots} | tr ':' ' '`
  ikp=1
  plots=''
  for drw in ${drwplts}; do
    if $(isInlist ${KindPlots} ${drw}); then
      case ${drw} in
        'map2Dsfc')
          plts=${pltdiffmap2Dsfc}
        ;;
        'map3D')
          plts=${pltdiffmap3D}
        ;;
        'shadcont2Dsfc')
          plts=${pltdiffshadcont2Dsfc}
        ;;
        'shadconthovmsfc')
          plts=${pltdiffshadconthovmsfc}
        ;;
        'shadcont2Dzsec')
          plts=${pltdiffshadcont2Dzsec}
        ;;
        *)
          echo ${errmsg}
          echo "  "${main}": plot kind '"${drw}"' not ready !!"
          exit
      esac
      # Do we have plots for this kind?
      Lkplots=`expr length ${plts}'0'`
      if test ${Lkplots} -lt 2; then
        ferrmsg 1 ${main} "plot!kind!"${drw}"!without!variables"
      fi
      kplots=`list_add_label ${plts} ${drw} + beg`
      if test ${ikp} -eq 1; then
        plots=${kplots}
      else
        plots=${plots}':'${kplots}
      fi
      ikp=`expr ${ikp} + 1`
    fi
# End of coincident direct plots
  done

  echo ${fname}" listplots= "${plots}_${ikp}
}

function draw_plot(){
# Function to carry out the drawing of the plots
#   plot= [kplot]+[varn1]|[kind1]#[varn2]|[kind2]#....
#     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
#
#     ':' separated list of statitsics variable values are given as: [var]|[kind(including combo values)]
#     in figures with more than 1 variable, use '#' to separate them
#   odir= working directory
#   headf= head of the file
#   additiona= additional value when necessary
#     'diff': diff plots (to use `diffspecificvarplot' instead)
  fname='draw_plot'

  plot=$1
  odir=$2
  headf=$3
  scratch=$4
  additional=$5

  if test $# -eq 5; then
    if test ${additional} = 'diff'; then
      specvarplot=${diffspecificvarplot}
    else
      specvarplot=${specificvarplot}
    fi
  fi

  cd ${odir}

  kplot=`echo ${plot} | tr '+' ' ' | awk '{print $1}'`
  plotvals=`echo ${plot} | tr '+' ' ' | awk '{print $2}'`

  # Assuming that multiple variables will be separated by '#'
  nvars=`echo ${plotvals} | tr '#' ' ' | wc -w | awk '{print $1}'`

  # Model and experiment
  mod=`echo ${odir} | tr '/' '\n' | tail -n 2 | head -n 1`
  exp=`echo ${odir} | tr '/' '\n' | tail -n 1`

  figfiles=''
  varinfilen=''
  CFvars=''
  figvarparam=''
  varkinds=''
  vardims=''
  if test ${nvars} -ge 1; then
    figname=${kplot}_${headf}
    iv=1
    while test $iv -le ${nvars}; do
      varkind=`echo ${plotvals} | tr '#' '\n' | head -n ${iv} | tail -n 1`
      varn=`echo ${varkind} | tr '|' ' ' | awk '{print $1}'`
      kind=`echo ${varkind} | tr '|' ' ' | awk '{print $2}'`

      # Do we have processed the given variable?
      nfiles=`ls -1 ${odir}/${varn}_${headf}*.nc | wc -l | awk '{print $1}'`
      if test ${nfiles} -eq 0; then
        echo "  "${fname}": there are no files for variable '"${varn}"' skiping it !!"
        return
      fi

      Nopers=`echo ${kind} | tr '@' ' ' | wc -w | awk '{print $1}'`
      if test ${Nopers} -eq 1; then
        if $(isInlist ${varmeanname} ${kind}); then 
          vn=${varn}'mean'
        else
          vn=${varn}
        fi
        if test ${op} = 'turb'; then
          vn=${vn}'turb'
        fi          
      else
        oprs=`echo ${kind} | tr '@' ' '`
        vn=${varn}
        for op in ${oprs}; do
          if $(isInlist ${varmeanname} ${op}); then 
            vn=${vn}'mean'
          fi
          if test ${op} = 'turb'; then
            vn=${vn}'turb'
          fi          
        done
      fi
      filen=`stats_filename ${varn} ${kind} ${headf}`
      dims=`python ${pyHOME}/nc_var.py -o idims -f ${filen} | grep alldims |         \
        awk '{print $3}'`

      if test ${iv} -eq 1; then
        if test ${Nopers} -gt 1; then
          varkinds=${varn}'|('$(echo ${kind} | tr '@' '|')')'
        else
          varkinds=${varkind}
        fi
        CFvars=${varn}
        varinfilen=${vn}
        figfiles=${filen}
        vardims=${dims}
        if $(isInlist ${dims} 'time'); then
          tunits=`python ${pyHOME}/nc_var.py -o ivattrs -f ${filen} -v time |        \
            grep units | grep -v '#' | awk '{print $3}'`
        fi
      else
        if test ${Nopers} -gt 1; then
          varkinds=${varkinds}':'${varn}'|('$(echo ${kind} | tr '@' '|')')'
        else
          varkinds=${varkinds}':'${varkind}
        fi
        CFvars=${CFvars}':'${varn}
        varinfilen=${varinfilen}':'${vn}
        figfiles=${figfiles}':'${filen}
        vardims=${vardims}':'${dims}
      fi

      figname=${figname}'_'${vn}'-'$(echo ${kind} | tr '@' '-')

      # Getting plot variable configuration (using variable name inside the file?)
      specific=false
#      if $(isInlist ${specvarplot} ${varn}); then
      if test $(index_string ${specvarplot} ${varn}) -ne -1; then
        vals=`list_filter ${specvarplot} ${varn}`
        kvals=`list_filter ${vals} ${kind}`
        Lkvals=`expr length ${kvals}'0'`
        if test ${Lkvals} -gt 1; then
          fkvals=`list_filter ${kvals} ${kplot}`
          Lfkvals=`expr length ${fkvals}'0'`
          if test ${Lfkvals} -gt 1; then
            specific=true
            minval=`echo ${fkvals} | tr '|' ' ' | awk '{print $4}'`
            maxval=`echo ${fkvals} | tr '|' ' ' | awk '{print $5}'`
            colorbar=`echo ${fkvals} | tr '|' ' ' | awk '{print $6}'`
            cntformat=`echo ${fkvals} | tr '|' ' ' | awk '{print $7}'`
            colorcnt=`echo ${fkvals} | tr '|' ' ' | awk '{print $8}' | tr ':' ' ' | awk '{print $1}'`
          fi
        fi
      fi
      if ! ${specific}; then
        allvals=`python ${pyHOME}/drawing.py -o variable_values -S ${varn} |         \
          grep all_values | awk '{print $3}'`
        pyn=$?
        ferrmsg ${pyn} ${fname} "python!'variable_values'!'${varn}'!failed#"
        minval=`echo ${allvals} | tr ',' ' ' | awk '{print $3}'`
        maxval=`echo ${allvals} | tr ',' ' ' | awk '{print $4}'`
        colorbar=`echo ${allvals} | tr ',' ' ' | awk '{print $7}'`
        cntformat='%g'
        colorcnt='black'
      fi

      # Graphical parameters for each variable in the plot
      if test ${iv} -eq 1; then
        figvarparam=${minval}'|'${maxval}'|'${colorbar}'|'${cntformat}'|'${colorcnt}
      else
        figvarparam=${figvarparam}':'${minval}'|'${maxval}'|'${colorbar}'|'
        figvarparam=${figvarparam}${cntformat}'|'${colorcnt}
      fi

      iv=`expr ${iv} + 1`
    # End of number of variables of the plot
    done
  fi
  figname=${figname}.${kindfig}

  if ${scratch}; then rm ${figname} >& /dev/null; fi
  if test ! -f ${figname}; then
    case ${kplot} in
      'diffmap2Dsfc')
        echo "  "${fname}": kind of plot '"${kplot}"' not ready !!"
        exit
      ;;
      'diffmap2Dz')
        echo "  "${fname}": kind of plot '"${kplot}"' not ready !!"
        exit
      ;;
      'map2Dsfc')
        echo "  "${fname}": kind of plot '"${kplot}"' not ready !!"
        exit
      ;;
      'map3D')
        echo "  "${fname}": kind of plot '"${kplot}"' not ready !!"
        exit
      ;;
      'shadcont2Dsfc')
        figtit=${mod}'|'${exp}'|'$(echo ${varkinds} | tr ':' '|')
        shdstdn=`echo ${CFvars} | tr ':' ' ' | awk '{print $1}'`
        cntstdn=`echo ${CFvars} | tr ':' ' ' | awk '{print $2}'`
        figfs=`echo ${figfiles} | tr ':' ','`
        srange=`echo ${figvarparam} | tr ':' ' ' | awk '{print $1}' | tr '|' ' ' |   \
          awk '{print $1","$2}'`
        cbar=`echo ${figvarparam} | tr ':' ' ' | awk '{print $1}' | tr '|' ' ' |     \
          awk '{print $3}'`
        crange=`echo ${figvarparam} | tr ':' ' ' | awk '{print $2}' | tr '|' ' ' |   \
          awk '{print $1","$2}'`
        cline=`echo ${figvarparam} | tr ':' ' ' | awk '{print $2}' | tr '|' ' ' |    \
          awk '{print $5}'`
        cfmt=`echo ${figvarparam} | tr ':' ' ' | awk '{print $2}' | tr '|' ' ' |     \
          awk '{print $4}'`

        graphvals=$(echo ${CFvars} | tr ':' ',')
        graphvals=${graphvals}':lon|-1,lat|-1:lon|-1,lat|-1:lon:lat:'${cbar}':fixc,'
        graphvals=${graphvals}${cline}':'${cfmt}':'${srange}':'${crange}',9:'
        graphvals=${graphvals}${figtit}':'${kindfig}':None:'${mapval}
 
        pyout=`python ${pyHOME}/drawing.py -f ${figfs} -o draw_2D_shad_cont          \
          -S ${graphvals} -v $(echo ${varinfilen} | tr ':' ',')`
        pyn=$?
        Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
        ferrmsg ${pyn} ${fname} "python!'draw_2D_shad_cont'!failed#"${Spyout}
        mv 2Dfields_shadow-contour.${kindfig} ${figname}

        # keeping all figures
        plotins="python "${pyHOME}"/drawing.py -f "${figfs}" -o draw_2D_shad_cont "
        plotins=${plotins}"-S '"${graphvals}"' -v "$(echo ${varinfilen} | tr ':' ',')
        echo " " >> ${odir}/all_figures.inf
        echo "#"$(echo $f{igtit} | tr '|' ' ') >> ${odir}/all_figures.inf
        echo ${plotins} >> ${odir}/all_figures.inf
      ;;
      'shadconthovmsfc')
        shdstdn=`echo ${CFvars} | tr ':' ' ' | awk '{print $1}'`
        cntstdn=`echo ${CFvars} | tr ':' ' ' | awk '{print $2}'`

        # It is assumed that if the space variable is 'lon': is desired a (lon, time) plot
        #   it it is 'lat': then (time, lat) plot
        dimspace='xxx' 
        if $( isInlist ${vardims} 'lon'); then
          spacedim='lon'
          figtit=${mod}'|'${exp}'|mean|longitudinal|evolution|of|'${shdstdn}'|&|'
          figtit=${figtit}${cntstdn}
          reverse='None'
          dims=${spacedim}'|-1,time|-1;'${spacedim}'|-1,time|-1;'${spacedim}';time'
        else
          spacedim='lat'
          figtit=${mod}'|'${exp}'|mean|meridional|evolution|of|'${shdstdn}'|&|'
          figtit=${figtit}${cntstdn}
          reverse='None'
          dims=${spacedim}'|-1,time|-1;'${spacedim}'|-1,time|-1;time;'${spacedim}
        fi

        figfs=`echo ${figfiles} | tr ':' ','`
        srange=`echo ${figvarparam} | tr ':' ' ' | awk '{print $1}' | tr '|' ' ' |   \
          awk '{print $1","$2}'`
        cbar=`echo ${figvarparam} | tr ':' ' ' | awk '{print $1}' | tr '|' ' ' |     \
          awk '{print $3}'`
        crange=`echo ${figvarparam} | tr ':' ' ' | awk '{print $2}' | tr '|' ' ' |   \
          awk '{print $1","$2}'`
        cline=`echo ${figvarparam} | tr ':' ' ' | awk '{print $2}' | tr '|' ' ' |    \
          awk '{print $5}'`
        cfmt=`echo ${figvarparam} | tr ':' ' ' | awk '{print $2}' | tr '|' ' ' |     \
          awk '{print $4}'`

        graphvals=$(echo ${CFvars} | tr ':' ',')
        graphvals=${shdstdn}','${cntstdn}';'${dims}';'
        graphvals=${graphvals}${cbar}';fixsigc,'${cline}';'${cfmt}';'${srange}';'
        graphvals=${graphvals}${crange}',9;'${figtit}';'${kindfig}';'${reverse}';'
        graphvals=${graphvals}'time|'${tunits}'|'${timekind}'|'${timefmt}'|'
        graphvals=${graphvals}${timelabel}

        python ${pyHOME}/drawing.py -f ${figfs} -o draw_2D_shad_cont_time            \
         -S ${graphvals} -v $(echo ${varinfilen} | tr ':' ',')
        mv 2Dfields_shadow-contour.${kindfig} ${figname}

        # keeping all figures
        plotins="python "${pyHOME}"/drawing.py -f "${figfs}" -o draw_2D_shad_cont_time"
        plotins=${plotins}" -S ${graphvals} -v "$(echo ${varinfilen} | tr ':' ',')
        echo " " >> ${odir}/all_figures.inf
        echo "#"$(echo $f{igtit} | tr '|' ' ') >> ${odir}/all_figures.inf
        echo ${plotins} >> ${odir}/all_figures.inf
      ;;
      'shadcont2Dzsec')
        shdstdn=`echo ${CFvars} | tr ':' ' ' | awk '{print $1}'`
        cntstdn=`echo ${CFvars} | tr ':' ' ' | awk '{print $2}'`
        # It is assumed that if the space variable is 'lon': is desired a (lon, pres) plot
        #   it it is 'lat': then (pres, lat) plot
        dimspace='xxx' 
        if $( isInlist ${vardims} 'lon'); then
          spacedim='lon'
          figtit=${mod}'|'${exp}'|mean|longitudinal|vertical|cross|section|of|'
          figtit=${figtit}${shdstdn}'|&|'${cntstdn}
        else
          spacedim='lat'
          figtit=${mod}'|'${exp}'|mean|meridional|vertical|cross|section|of|'
          figtit=${figtit}${shdstdn}'|&|'${cntstdn}
        fi
        reverse='flip@y'
        dims=${spacedim}'|-1,pres|-1:'${spacedim}'|-1,pres|-1:'${spacedim}':pres'

        figtit=${mod}'|'${exp}'|'$(echo ${varkinds} | tr ':' '|')
        figfs=`echo ${figfiles} | tr ':' ','`
        srange=`echo ${figvarparam} | tr ':' ' ' | awk '{print $1}' | tr '|' ' ' |   \
          awk '{print $1","$2}'`
        cbar=`echo ${figvarparam} | tr ':' ' ' | awk '{print $1}' | tr '|' ' ' |     \
          awk '{print $3}'`
        crange=`echo ${figvarparam} | tr ':' ' ' | awk '{print $2}' | tr '|' ' ' |   \
          awk '{print $1","$2}'`
        cline=`echo ${figvarparam} | tr ':' ' ' | awk '{print $2}' | tr '|' ' ' |    \
          awk '{print $5}'`
        cfmt=`echo ${figvarparam} | tr ':' ' ' | awk '{print $2}' | tr '|' ' ' |     \
          awk '{print $4}'`

        graphvals=$(echo ${CFvars} | tr ':' ',')
        graphvals=${graphvals}':'${dims}':'${cbar}':fixc,'
        graphvals=${graphvals}${cline}':'${cfmt}':'${srange}':'${crange}',9:'
        graphvals=${graphvals}${figtit}':'${kindfig}':'${reverse}':None'
 
        pyout=`python ${pyHOME}/drawing.py -f ${figfs} -o draw_2D_shad_cont          \
          -S ${graphvals} -v $(echo ${varinfilen} | tr ':' ',')`
        pyn=$?
        Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
        ferrmsg ${pyn} ${fname} "python!'draw_2D_shad_cont'!failed#"${Spyout}
        mv 2Dfields_shadow-contour.${kindfig} ${figname}

        # keeping all figures
        plotins="python "${pyHOME}"/drawing.py -f "${figfs}" -o draw_2D_shad_cont "
        plotins=${plotins}"-S '"${graphvals}"' -v "$(echo ${varinfilen} | tr ':' ',')
        echo " " >> ${odir}/all_figures.inf
        echo "#"$(echo $f{igtit} | tr '|' ' ') >> ${odir}/all_figures.inf
        echo ${plotins} >> ${odir}/all_figures.inf
      ;;
      *)
        echo ${errormsg}
        echo "  "${fname}": plot kind '"${kplot}"' not ready !!"
        exit
      ;;
    esac
  fi

#LLUIS

}

#######    #######
## MAIN
    #######
rootsh=`pwd`

# Uploading environment
uploadvars model_graphics.dat

if test ${scratch} = 'true'; then
  scratch=true
  echo ${warnmsg}
  echo "  "${main}": starting from the SCRATCH !!"
  echo "    10 seconds left!!"
  filescratch=true
  figscratch=true
  sleep 10
else
  scratch=false
  if test ${filescratch} = 'true'; then
    filescratch=true
    echo ${warnmsg}
    echo "  "${main}": files starting from the SCRATCH !!"
    echo "    5 seconds left!!"
    sleep 5
  else
    filescratch=false
  fi
  if test ${figscratch} = 'true'; then
    figscratch=true
    echo ${warnmsg}
    echo "  "${main}": figures starting from the SCRATCH !!"
    echo "    5 seconds left!!"
    sleep 5
  else
    figscratch=false
  fi
fi

if test ${addfiles} = 'true'; then
  addfiles=true
else
  addfiles=false
fi

if test ${addfigures} = 'true'; then
  addfigures=true
else
  addfigures=false
fi

if test ${debug} = 'true'; then
  dbg=true
else
  dbg=false
fi

timeval='tstep'
zval='null'
dimval='lon,lat,pressure,time'

compute=0
computetlon=false
plot=false
single2Dplot=false
lonZsec=false
differences=true

combosfile=${pyHOME}/diagnostics.inf

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

mods=`echo ${models} | tr ':' ' '`
varks=`echo ${varkinds} | tr ':' ' '`

# Models loop
##
for mod in ${mods}; do
  case ${mod} in
    'WRF')
      exps=`echo ${WRFexps} | tr ':' ' '`
      fheaders=`echo ${WRFheaders} | tr ':' ' '`
    ;;
    'LMDZ')
      exps=`echo ${LMDZexps} | tr ':' ' '`
      fheaders=`echo ${LMDZheaders} | tr ':' ' '`
    ;;
    'WRF_LMDZ')
      exps=`echo ${WRF_LMDZexps} | tr ':' ' '`
      fheaders=`echo ${WRF_LMDZheaders} | tr ':' ' '`
    ;;
    '*')
      echo ${errmsg}
      echo "  "${main}": model '"${mod}"' not ready!!"
      exit
    ;;
  esac  

  modinf=`${pyHOME}/nc_var.py -o model_characteristics -f None -S ${mod} | \
    grep singleline | awk '{print $2}'`
  ferrmsg $? $main "python!'model_characteristics'!failed"
  dnx=`echo ${modinf} | tr ';' '\n' | grep dimxn | tr '=' ' ' | awk '{print $2}'`
  dny=`echo ${modinf} | tr ';' '\n' | grep dimyn | tr '=' ' ' | awk '{print $2}'`
  dnz=`echo ${modinf} | tr ';' '\n' | grep dimzn | tr '=' ' ' | awk '{print $2}'`
  dnt=`echo ${modinf} | tr ';' '\n' | grep dimtn | tr '=' ' ' | awk '{print $2}'`
  vdnx=`echo ${modinf} | tr ';' '\n' | grep vardxn | tr '=' ' ' | awk '{print $2}'`
  vdny=`echo ${modinf} | tr ';' '\n' | grep vardyn | tr '=' ' ' | awk '{print $2}'`
  vdnz=`echo ${modinf} | tr ';' '\n' | grep vardzn | tr '=' ' ' | awk '{print $2}'`
  vdnt=`echo ${modinf} | tr ';' '\n' | grep vardtn | tr '=' ' ' | awk '{print $2}'`
  echo ${mod}
  echo "  dims: "${dnx}", "${dny}", "${dnz}", "${dnt}
  echo "  var dims: "${vdnx}", "${vdny}", "${vdnz}", "${vdnt}
  moddims=${dnx}','${dny}','${dnz}','${dnt}
  modvdims=${vdnx}','${vdny}','${vdnz}','${vdnt}

# Experiments loop
##
  for exp in ${exps}; do
    echo "  "${exp}"..."
    iwdir=${ifold}/${mod}/${exp}

    # Does input folder exist?
    if test ! -d ${iwdir}; then
      echo ${errmsg}
      echo "  "${main}": folder '"${iwdir}"' does not exist !!"
      exit
    fi

    owdir=${ofold}/${mod}/${exp}
    mkdir -p ${owdir}

    # Need to pass to analyze all the data?
    if ${filescratch}; then 
      rm ${owdir}/varcompute.inf >& /dev/null
      rm ${owdir}/all_computevars.inf >& /dev/null
      rm ${owdir}/all_statsvars.inf >& /dev/null

      echo "## Computation of variables " > ${owdir}/all_computevars.inf
      echo "## Computation of statistics " > ${owdir}/all_statsvars.inf
    fi

    if ${addfiles}; then
      rm ${owdir}/varcompute.inf >& /dev/null
    fi

    if test ! -f ${owdir}/varcompute.inf; then

      # Does input folder has header files?
      cd ${iwdir}
      files=''
      testfiles=''
      ih=1
      for fh in ${fheaders}; do
        if ${filescratch}; then rm ${owdir}/*_${fh}*.nc >& /dev/null; fi

        files1h=`ls -1 ${fh}* | tr '\n' '@'`
        Lfiles1h=`expr length ${files1h}'0'`
        if test ${Lfiles1h} -lt 2; then
          ferrmsg 1 $main "folder!:!"${iwdir}"!does!not!contain!files!"${fh}"*"
        fi
        testfiles1h=`ls -1 ${fh}* | head -n 1`
        if test ${ih} -eq 1; then
          files=${fh}'|'${files1h}
          testfiles=${fh}'|'${testfiles1h}
        else
          files=${files}','${fh}'|'${files1h}
          testfiles=${testfiles}'@'${fh}'|'${testfiles1h}
        fi
        ih=`expr ${ih} + 1`
      done
      testfs=`echo ${testfiles} | tr '@' ' '`
      cd ${rootsh}

      varcompt=`compvars_listconstruct ${varkinds} $(echo ${testfiles} | tr '@' ':') \
       | grep varcompute | awk '{print $3}'`

      varcompute=`echo ${varcompt} | tr '_' ' ' | awk '{print $1}'`
      itotv=`echo ${varcompt} | tr '_' ' ' | awk '{print $2}'`
      combo=`echo ${varcompt} | tr '_' ' ' | awk '{print $3}'`

      # Outwritting the varcompute to avoid next time (if it is not filescratch!)
      cat << EOF > ${owdir}/varcompute.inf
files: ${files}
varcompute: ${varcompute}
itotv: ${itotv}
testfs: ${testfiles}
combo: ${combo}
EOF
    else
      echo $warnmsg
      echo "  "${main}": getting already data information from the experiment!"
      files=`cat ${owdir}/varcompute.inf | grep files | awk '{print $2}'`
      varcompute=`cat ${owdir}/varcompute.inf | grep varcompute | awk '{print $2}'`
      itotv=`cat ${owdir}/varcompute.inf | grep itotv | awk '{print $2}'`
      testfiles=`cat ${owdir}/varcompute.inf | grep testfs | awk '{print $2}'`
      combo=`cat ${owdir}/varcompute.inf | grep 'combo:' | awk '{print $2}'`
# End of avoiding to repeat all the experiment search
    fi

    echo "  For experiment '"${exp}"' is required to compute: "${itotv}" variables"
# Computing files for each variable
## 
    echo "      Computing all variables ..."

    compute_vars ${iwdir} ${owdir} ${files} ${moddims} ${modvdims} ${varcompute}     \
      ${filescratch}

    # Computing combos
    if test ${combo}; then
      echo "      Computing combos: "${varcombo}
      compute_combos ${iwdir} ${owdir} ${files} ${testfiles} ${moddims} ${modvdims}  \
        ${varcombo} ${filescratch}
    fi

    echo "  "${main}": "${ic}" variables has been computed"
    echo "  "${main}": "${isc}" statistics has been computed"

# Direct plotting
##
    echo "  "${main}": Plotting direct figures ..."
    if ${figscratch}; then rm directplotsdraw.inf; fi
    if ${addfigures}; then rm directplotsdraw.inf; fi

    if test ! -f directplotsdraw.inf; then
      listplots=`createlist_plots ${drawplots} ${directplots} | grep listplots | awk '{print $3}'`

      plots=`echo ${listplots} | tr '_' ' ' | awk '{print $1}'`
      ikp=`echo ${listplots} | tr '_' ' ' | awk '{print $2}'`

      cat << EOF > directplotsdraw.inf
Nkinds= ${ikp}
plots= ${plots}
EOF
    else
      echo ${warnmsg}
      echo "  "${main}": retrieving direct plots from existing file 'directplotsdraw.inf' !!"
      Nkinds=`cat directplotsdraw.inf | grep Nkinds | awk '{print $2}'`
      plots=`cat directplotsdraw.inf | grep plots | awk '{print $2}'`
    fi

    if ${figscratch}; then rm ${owdir}/all_figures.inf; fi
    if test ! -f ${owdir}/all_figures.inf; then
      echo "###### List of all figures" > ${owdir}/all_figures.inf
    fi

    pts=`echo ${plots} | tr ':' ' '`
    idp=1
    for pt in ${pts}; do
      echo "        "${pt}
      draw_plot ${pt} ${owdir} ${fileh} ${figscratch}
      idp=`expr ${idp} + 1`
    done

    echo "  "${main}": "${idp}" direct plots have been computed"
    cd ${rootsh}

# end of experiments
  done

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

  diffexperiments=`list_pairs $(echo ${exps} | tr ' ' ':') | grep -v newlists | awk '{print $3}'`
  diffexps=`echo ${diffexperiments} | tr ':' ' '`

  for diffexp in ${diffexps}; do
    echo "  "${diffexp}
    difexp1=`echo ${diffexp} | tr '|' ' ' | awk '{print $1}'`
    difexp2=`echo ${diffexp} | tr '|' ' ' | awk '{print $2}'`

    vdiffs=`echo ${vardiff} | tr ':' ' '`
    ivd=1
    echo "      Computing diff statistics ..."
    for vdiff in ${vdiffs}; do
      echo ${vdiff}" ..."
      if test $(expr index ${vdiff} '|') -ne -1; then
        vn=`echo ${vdiff} | tr '|' ' ' | awk '{print $1}'`
        ops=`echo ${vdiff} | tr '|' ' ' | awk '{print $2}'`
      else
        vn=`echo ${vdiff} | tr '|' ' ' | awk '{print $1}'`
        ops='None'
      fi
      #Modify in case of combo
      if test $(expr index ${ops} '@') -eq 0; then
        fileh=`echo ${varcompute} | tr ';' '\n' | grep ${vdiff} | tr '|' ' ' |       \
          awk '{print $3}'`
      else
        set -x
        Lcombo=`expr length ${vn}';'${ops}'|combo|'`
        Lvarcompute=`expr length ${varcompute}`
        indcombo=`index_string ${varcompute} ${vn}';'${ops}'|combo|'`
        Lindcombo=`expr ${Lvarcompute} - ${indcombo}`

        fileh=`echo ${varcompute:${indcombo}:${Lindcombo}} | tr ';' '\n' |           \
          head -n 2  | tail -n 1 | tr '|' ' ' | awk '{print $3}'`
        if test $(index_string ${ops} pinterp) -ne -1; then fileh=${fileh}'p_pinterp'; fi
      fi

      owdir=${ofold}/${mod}/${difexp1}-${difexp2}
      if test ${ivd} -eq 1; then
        mkdir -p ${owdir}
      fi

      filen=${vn}_${fileh}.nc

      if ${filescratch}; then rm ${owdir}/${filen} >& /dev/null; fi
      if test ! -f ${owdir}/${filen}; then
        set -x
        dfile1=${ofold}/${mod}/${difexp1}/${filen}
        dfile2=${ofold}/${mod}/${difexp2}/${filen}

        dfiles='add|'${dfile1}'|'${vn}',sub|'${dfile2}'|'${vn}
        dims='time|pres|lat|lon@'${dfiles}
        Lfileh=`expr length ${fileh}`
        Lfileh1=`expr ${Lfileh} - 1`
        echo "python ${pyHOME}/nc_var.py -S ${dims} -o compute_opersvarsfiles -v ${vn}"
        pyout=`python ${pyHOME}/nc_var.py -S ${dims} -o compute_opersvarsfiles -v ${vn}`
        pyn=$?
        Spyout=`echo ${pyout} | tr '\n' '#' | tr ' ' '!'`
        ferrmsg ${pyn} ${fname} "python!'difference'!failed"${Spyout}

        mv opersvarsfiles_${vn}.nc ${owdir}/${filen}
      fi

      cfdims='lon:lat:pres:time'
      cfvdims='lon:lat:pres:time'
      cdiffvn=`echo ${filen} | tr '_' ' ' | awk '{print $1}'`
      cdiffh=`echo ${filen} | tr '_' ' ' | awk '{print $2}' | tr '.' ' ' |           \
        awk '{print $1}'`

      cdiffvar=${cdiffvn}'|'${ops}'|'${cdiffh}'|None|None'

      if test $(expr index ${ops} '@') -eq 0; then
        compute_statistics ${owdir} ${filen} ${owdir} ${cdiffvar} ${cfdims}          \
          ${cfvdims} ${filescratch}
      else
        # Computing combos
        echo "      Computing combos: "${ops}
        combs=`echo ${ops} | tr '@' ' '`
        ifile=${cdiffvn}_${fileh}
        for comb in ${combs}; do
          combcvar=${cdiffvn}'|'${comb}'|'${cdiffh}'|None|None'
          compute_statistics ${owdir} ${ifile}.nc ${owdir} ${combcvar} ${cfdims}     \
            ${cfvdims} ${filescratch} ${ifile}.nc
          ifile=${ifile}_${comb}

          # Adding 'surnames' to the variable's name
          if test ${comb} = 'turb'; then
            cdiffvn=${cdiffvn}'turb'
          fi
          if $(isInlist ${varmeanname} ${comb}); then
            cdiffvn=${cdiffvn}'mean'
          fi

        # End of combo 
        done
      fi
      ivd=`expr ${ivd} + 1`
    # end of variable diff
    done

# Diff plotting
##
    echo "  "${main}": Plotting diff figures ..."
    if ${figscratch}; then rm diffplotsdraw.inf; fi
    if ${addfigures}; then rm diffplotsdraw.inf; fi

    if test ! -f diffplotsdraw.inf; then
      listplots=`createlist_diffplots ${drawdiffplots} ${directplots} |              \
        grep listplots | awk '{print $3}'`

      plots=`echo ${listplots} | tr '_' ' ' | awk '{print $1}'`
      ikp=`echo ${listplots} | tr '_' ' ' | awk '{print $2}'`

      cat << EOF > diffplotsdraw.inf
Nkinds= ${ikp}
plots= ${plots}
EOF
    else
      echo ${warnmsg}
      echo "  "${main}": retrieving direct plots from existing file 'diffplotsdraw.inf' !!"
      Nkinds=`cat diffplotsdraw.inf | grep Nkinds | awk '{print $2}'`
      plots=`cat diffplotsdraw.inf | grep plots | awk '{print $2}'`
    fi

    pts=`echo ${plots} | tr ':' ' '`
    idp=1
    for pt in ${pts}; do
      echo "        "${pt}
      draw_plot ${pt} ${owdir} ${fileh} ${figscratch} 'diff'
      idp=`expr ${idp} + 1`
    done

    echo "  "${main}": "${idp}" diff plots have been computed"
    cd ${rootsh}

  # end experiment pairs
  done
# LLUIS

  exit

# end of models
done

exit

istep=0

if test ${compute} -eq 1; then
  for exp in ${exps}; do
    echo "exp: "${exp}
    if test ${exp} == 'AR40'; then
      var2Dtmeans=${var2Dtmeans_AR40}
    else
      var2Dtmeans=${var2Dtmeans_NPv31}
    fi
    for vark in ${varks}; do
      echo "  vark: "${vark}

      if test ${vark} = 'mean'; then
        fileorig=${ifoldmean}/${exp}/${filenmean}
        if test ! -f ${fileorig}; then
          fileworig=${filensfc}
          dvals='T:Time,Z:bottom_top,Y:south_north,X:west_east'
          dimnames='T:Time,Z:bottom_top,Y:south_north,X:west_east'
          python ${pyHOME}/vertical_interpolation.py                            \
            -f ${ofold}/${exp}/${filensfc} -o WRFp -i ${plevels} -k 'lin'            \
            -v WRFt,U,V,WRFrh,WRFght -d ${dimnames} -D T:Times,Z:ZNU,Y:XLAT,X:XLONG
          if test $? -ne 0; then
            echo ${errormsg}
            echo "  python failed!"
            echo "  python ${pyHOME}/vertical_interpolation.py"                            \
            "-f ${ofold}/${exp}/${filensfc} -o WRFp -i ${plevels} -k 'lin'"            \
            "-v WRFt,U,V,WRFrh,WRFght -d ${dimnames} -D T:Times,Z:ZNU,Y:XLAT,X:XLONG"
            rm ${ofold}/${exp}/vertical_interpolation_WRFp.nc >& /dev/null
            exit
          fi
          mv vertical_interpolation_WRFp.nc ${ofold}/${exp}
        fi
        vars=`echo ${mean_variables} | tr ':' ' '`
        ofile='mean_variables'
      elif test ${vark} = 'sfc'; then
        fileorig=${ifoldsfc}/${exp}/${filensfc}
        if test ! -f ${fileorig}; then
          fileworig=${filensfc}
          cp ${ifoldsfc}/${exp}/${fileworig} ${ifoldsfc}/${exp}/${filensfc}
          python ${pyHOME}/nc_var.py -o WRF_CFxtime_creation -f ${ifoldsfc}/${exp}/${filensfc} \
            -S 19491201000000,hours -v time
          python ${pyHOME}/nc_var.py -o WRF_CFlonlat_creation -f ${ifoldsfc}/${exp}/${filensfc}\
            -S longitude,latitude -v XLONG,XLAT
        fi
        if test ${exp} = 'AR40'; then
          vars=`echo ${sfc_variables_AR40} | tr ':' ' '`
        else
          vars=`echo ${sfc_variables_NPv31} | tr ':' ' '`
        fi
        ofile='sfc_variables'
      elif test ${vark} = 'diag'; then
        fileorig=${ifolddiag}/${exp}/${filendiag}
        vars=`echo ${diag_variables} | tr ':' ' '`
        if test ! -f ${fileorig}; then
          values='Time@XTIME,bottom_top@ZNU,south_north@XLAT,west_east@XLONG'
          ivar=0
          varcombos=''
          for var in ${vars}; do
            echo "    var:"${var}
            varn=`echo ${var} | tr '_' ' ' | awk '{print $2}'`
            combo=`python ${pyHOME}/diagnostics.py -f ${combosfile} -d variable_combo -v ${varn}\
              | grep COMBO | awk '{print $2}'`
            if test ${combo} = 'ERROR'; then
              echo ${errormsg}
              echo "  No variable '"${varn}"' in '"${combosfile}"' !!"
              exit
            fi
            if test ${ivar} -eq 0; then
              varcombos=${varn}'|'${combo}
            else
              varcombos=${varcombos}','${varn}'|'${combo}
            fi
            ivar=`expr ${ivar} + 1`
          done
          python ${pyHOME}/diagnostics.py -d ${values} -f ${ifoldsfc}/${exp}/${filensfc} -v \
            ${varcombos}
          if test $? -ne 0; then
            echo ${errormsg}
            echo "  python failed!!"
            echo python ${pyHOME}/diagnostics.py -d ${values} -f ${ifoldsfc}/${exp}/${filensfc}\
              -v ${varcombos}
            exit
          fi
          python ${pyHOME}/nc_var.py -o WRF_CFxtime_creation -f diagnostics.nc                 \
            -S 19491201000000,hours -v time
           python ${pyHOME}/nc_var.py -o WRF_CFlonlat_creation -f diagnostics.nc               \
            -S longitude,latitude -v XLONG,XLAT
          mv diagnostics.nc ${ifolddiag}/${exp}
        fi
        ofile='diag_variables'
      elif test ${vark} = 'z'; then
        fileorig=${ifoldmean}/${exp}/${filenmean}
        vars=`echo ${z_variables} | tr ':' ' '`
        ofile='z_variables'
      fi

# Averaging
##
      echo "  averaging..."
      ivar=0
      for varn in ${vars}; do
        file=${fileorig}
        echo "    var: "${varn}
        var=`echo ${varn} | tr '|' ' ' | awk '{print $2}'`
        varval=`python ${pyHOME}/drawing.py -o variable_values -S ${var} | grep all_values | awk '{print $3}'`
        varname=`echo ${varval} | tr ',' ' ' | awk '{print $1}'`

        echo "      std name: "${varname}
        varhead=`echo ${varname} | tr '_' ' ' | awk '{print $1}'`
        vartail=`echo ${varname} | tr '_' ' ' | awk '{print $2}'`

        if test ${vark} = 'mean'; then
          vals='time:-1|z:-1|y:-1|x:-1,time:x,mean,pressure:lat'
          python ${pyHOME}/nc_var.py -o file_oper_alongdims -S ${vals} -f ${file} -v ${var}
          pyexec=$?
          if test ${pyexec} -ne 0; then
            echo ${errormsg}
            echo "  "${main}": python fai1ls!"
            echo "python ${pyHOME}/nc_var.py -o file_oper_alongdims -S ${vals} -f ${file} -v ${var}"
            exit
          else
            oper=`echo ${vals} | tr ',' ' ' | awk '{print $3}'`
            mv file_oper_alongdims_mean.nc ${ofold}/${exp}/${varname}${oper}.nc
          fi
        elif test ${vark} = 'z'; then
          vals='time:-1|z:-1|y:-1|x:-1,time:x,mean,pressure:lat'
          echo "NOT FINISHED!!!!"
          exit
        else
          vals='Time:-1|south_north:-1|west_east:-1,west_east,mean,time:XLAT'
          if test ${var} = 'ACRAINTOT'; then
            files='add|'${file}'|RAINC,add|'${file}'|RAINNC'
            python ${pyHOME}/nc_var.py -S 'time|XLAT|XLONG@'${files} -o compute_opersvarsfiles \
              -v ${var}
            mv opersvarsfiles_${var}.nc ${ofold}/${exp}
            file=${ofold}/${exp}/opersvarsfiles_${var}.nc
          elif test ${var} = 'RAINTOT'; then
            dims='Time@time,bottom_top@ZNU,south_north@XLAT,west_east@XLONG'
            python ${pyHOME}/diagnostics.py -d ${dims} -v 'RAINTOT|RAINC@RAINNC@time' -f ${file}
            mv diagnostics.nc ${ofold}/${exp}/diagnostics_${varname}.nc
            file=${ofold}/${exp}/diagnostics_${varname}.nc
            var='pr'
          elif test ${var} = 'cllmh'; then
            dims='Time@time,bottom_top@ZNU,south_north@XLAT,west_east@XLONG'
            clouds='cll:clm:clh'
            cls=`echo ${clouds} | tr ':' ' '`
            for cl in ${cls}; do
              file=${ofold}/${exp}/diagnostics.nc
              var=${cl}
              echo "    var: "${var}
              python ${pyHOME}/nc_var.py -o file_oper_alongdims -S ${vals} -f ${file} -v ${cl}
              pyexec=$?
              if test ${pyexec} -ne 0; then
                echo ${errormsg}
                echo "  "${main}": python fails!"
                echo python ${pyHOME}/nc_var.py -o file_oper_alongdims -S ${vals} -f ${file} -v ${cl}
                exit
              fi
              mv file_oper_alongdims_mean.nc ${ofold}/${exp}/${cl}${oper}.nc
              vals='Time:-1|south_north:-1|west_east:-1,Time,mean,XLONG:XLAT'
              dims='Time@time,bottom_top@ZNU,south_north@XLAT,west_east@XLONG'
              python ${pyHOME}/nc_var.py -o file_oper_alongdims -S ${vals} -f ${file} -v ${cl}
              pyexec=$?
              if test ${pyexec} -ne 0; then
                echo ${errormsg}
                echo "  "${main}": python fails!"
                exit
              else
                oper=`echo ${vals} | tr ',' ' ' | awk '{print $3}'`
                mv file_oper_alongdims_mean.nc ${ofold}/${exp}/${cl}t${oper}.nc
              fi
            done
#            break
          elif test ${var} = 'WRFbils' && test ! -f ${ofold}/${exp}/diagnostics_${varname}.nc; then
            dims='Time@time,south_north@XLAT,west_east@XLONG'
            python ${pyHOME}/diagnostics.py -d ${dims} -v 'WRFbils|HFX@LH' -f ${file}
            mv diagnostics.nc ${ofold}/${exp}/diagnostics_${varname}.nc
            file=${ofold}/${exp}/diagnostics_${varname}.nc
            var='bils'
          elif test ${var} = 'WRFprw' && test ! -f ${ofold}/${exp}/diagnostics_${varname}.nc; then
            dims='Time@time,bottom_top@ZNU,south_north@XLAT,west_east@XLONG'
            python ${pyHOME}/diagnostics.py -d ${dims} -v 'WRFprw|WRFdens@QVAPOR' -f ${file}
            mv diagnostics.nc ${ofold}/${exp}/diagnostics_${varname}.nc
            file=${ofold}/${exp}/diagnostics_${varname}.nc
            var='prw'
          elif test ${var} = 'WRFrhs' && test ! -f ${ofold}/${exp}/diagnostics_${varname}.nc; then
            dims='Time@time,bottom_top@ZNU,south_north@XLAT,west_east@XLONG'
            python ${pyHOME}/diagnostics.py -d ${dims} -v 'WRFrhs|PSFC@T2@Q2' -f ${file}
            mv diagnostics.nc ${ofold}/${exp}/diagnostics_${varname}.nc
            file=${ofold}/${exp}/diagnostics_${varname}.nc
            var='huss'
          elif test ${var} = 'WRFprc'; then
            vals='Time:-1|south_north:-1|west_east:-1,west_east,mean,time:XLAT'
            dims='Time@time,bottom_top@ZNU,south_north@XLAT,west_east@XLONG'
            var='prc'
          elif test ${var} = 'WRFprls'; then
            vals='Time:-1|south_north:-1|west_east:-1,west_east,mean,time:XLAT'
            dims='Time@time,bottom_top@ZNU,south_north@XLAT,west_east@XLONG'
            var='prls'
          fi
          if test ! -f ${ofold}/${exp}/${varname}${oper}.nc; then
            python ${pyHOME}/nc_var.py -o file_oper_alongdims -S ${vals} -f ${file} -v ${var}
            pyexec=$?
            if test ${pyexec} -ne 0; then
              echo ${errormsg}
              echo "  "${main}": python fails!"
              exit
            else
              mv file_oper_alongdims_mean.nc ${ofold}/${exp}/${varname}${oper}.nc
            fi
          fi
        fi

# Time means
        if $(isin_list ${var2Dtmeans} ${var}); then 
          echo "Computing time means!!"
          if test ${var} = 'cllmh'; then
            clouds='cll:clm:clh'
            clds=`echo ${clouds} | tr ':' ' '`
            for cld in ${clds}; do
              vals='Time:-1|south_north:-1|west_east:-1,Time,mean,XLONG:XLAT'
              dims='Time@time,bottom_top@ZNU,south_north@XLAT,west_east@XLONG'
              python ${pyHOME}/nc_var.py -o file_oper_alongdims -S ${vals} -f ${file} -v ${cld}
              pyexec=$?
              if test ${pyexec} -ne 0; then
                echo ${errormsg}
                echo "  "${main}": python fails!"
                exit
              else
                oper=`echo ${vals} | tr ',' ' ' | awk '{print $3}'`
                mv file_oper_alongdims_mean.nc ${ofold}/${exp}/${cld}t${oper}.nc
              fi
            done
          else
            vals='Time:-1|south_north:-1|west_east:-1,Time,mean,XLONG:XLAT'
            dims='Time@time,bottom_top@ZNU,south_north@XLAT,west_east@XLONG'
            python ${pyHOME}/nc_var.py -o file_oper_alongdims -S ${vals} -f ${file} -v ${var}
            pyexec=$?
            if test ${pyexec} -ne 0; then
              echo ${errormsg}
              echo "  "${main}": python fails!"
              exit
            else
              oper=`echo ${vals} | tr ',' ' ' | awk '{print $3}'`
              mv file_oper_alongdims_mean.nc ${ofold}/${exp}/${varname}t${oper}.nc
            fi
          fi
        fi
        if test ${var} = 'cllmh'; then exit; fi
        ivar=`expr ${ivar} + 1`
#        exit
# end of vars
      done

# end of kind vars
    done

  done
# end of compute
fi
# Longitudinal means of tmeans
##
if ${computetlon}; then
  echo "Computing Longitudinal means of temporal means..."
  for exp in ${exps}; do
    echo "  exp:"${exp}
    for tmean in ${ofold}/${exp}/*tmean.nc; do
      fname=`basename ${tmean}`
      var=`ncdump -h ${tmean} | grep float | grep mean | tr '(' ' ' | awk '{print $2}'`
      vals='south_north:-1|west_east:-1,west_east,mean,XLAT'
      dims='bottom_top@bottom_top,south_north@south_north,west_east@west_east'
      if test ! -f ${ofold}/${exp}/${var}t-lonmean.nc; then
        python ${pyHOME}/nc_var.py -o file_oper_alongdims -S ${vals} -f ${tmean} -v ${var}
        pyexec=$?
        if test ${pyexec} -ne 0; then
          echo ${errormsg}
          echo "  "${main}": python fails!"
          exit
        else
          mv file_oper_alongdims_mean.nc ${ofold}/${exp}/${var}t-lonmean.nc
        fi
#      exit
      fi
    done

  done
fi

if ${plot}; then
# Plots
##

for exp in ${exps}; do
  echo "exp: "${exp}
  for vark in ${varks}; do
    echo "  vark: "${vark}
    if test ${vark} = 'mean'; then
      gcoups=`echo ${graph_couples_mean} | tr ':' ' '`
    elif test ${vark} = 'sfc'; then
      if test ${exp} = 'AR40'; then
        gcoups=`echo ${graph_couples_sfc_AR40} | tr ':' ' '`
      else
        gcoups=`echo ${graph_couples_sfc_NPv31} | tr ':' ' '`
      fi
    elif test ${vark} = 'diag'; then
      gcoups=`echo ${graph_couples_diag} | tr ':' ' '`
    fi

    for gpair in ${gcoups}; do
      shdvar=`echo ${gpair} | tr '@' ' ' | awk '{print $1}'`
      cntvar=`echo ${gpair} | tr '@' ' ' | awk '{print $2}'`
      echo "  couple: "${shdvar}'-'${cntvar}
      shdvals=`python ${pyHOME}/drawing.py -o variable_values -S ${shdvar} | grep all_values | awk '{print $3}'`
      cntvals=`python ${pyHOME}/drawing.py -o variable_values -S ${cntvar} | grep all_values | awk '{print $3}'`
      files=${ofold}'/'${exp}'/'${shdvar}'mean.nc,'${ofold}'/'${exp}'/'${cntvar}'mean.nc'

      Lshdvals=`expr length ${shdvals}0`
      if test ${Lshdvals} -lt 2; then
        echo ${errormsg}
        echo "  Error in drawing_tools.py 'variables_values'!!"
        echo "    variable '"${shdvar}"' NOT found!"
        exit
      fi

      Lcntvals=`expr length ${cntvals}0`
      if test ${Lcntvals} -lt 2; then
        echo ${errormsg}
        echo "  Error in drawing_tools.py 'variables_values'!!"
        echo "    variable '"${cntvar}"' NOT found!"
        exit
      fi

      shdstdn=`echo ${shdvals} | tr ',' ' ' | awk '{print $1}'`
      shdcbar=`echo ${shdvals} | tr ',' ' ' | awk '{print $7}'`
      shdmin=`echo ${shdvals} | tr ',' ' ' | awk '{print $3}'`
      shdmax=`echo ${shdvals} | tr ',' ' ' | awk '{print $4}'`

      cntstdn=`echo ${cntvals} | tr ',' ' ' | awk '{print $1}'`
      cntmin=`echo ${cntvals} | tr ',' ' ' | awk '{print $3}'`
      cntmax=`echo ${cntvals} | tr ',' ' ' | awk '{print $4}'`

      if test ${shdstdn} = 'ERROR' || test ${cntstdn} = 'ERROR'; then
        echo ${errmsg}
        echo "  "${main}": wrong variable names!!!"
        echo "    shdvals: "${shdvals}
        echo "    cntvals: "${cntvals}
        echo "    shdvar: "${shdvar}" cntvar: "${cntvar}
        exit
      fi

      cntcolor='black'
      cntfmt='%g'

      if test ${gpair} = 'va@ua'; then
        shdmin=`echo ${shdmin} | awk '{print $1/4.}'`
        shdmax=`echo ${shdmax} | awk '{print $1/4.}'`
        cntmin=`echo ${cntmin} | awk '{print $1/3.}'`
        cntmax=`echo ${cntmax} | awk '{print $1 + 20.}'`
      elif test ${gpair} = 'tas@ps'; then
        shdmin=`echo ${shdmin} | awk '{print $1 + 50.}'`
        shdmax=`echo ${shdmax} | awk '{print $1 - 15.}'`
        cntmin=`echo ${cntmin} | awk '{print $1 + 15000.}'`
        cntmax=`echo ${cntmax} | awk '{print $1 - 2000.}'`
      elif test ${gpair} = 'uas@vas'; then
        cntmin=`echo ${cntmin} | awk '{print $1 * 0.1}'`
        cntmax=`echo ${cntmax} | awk '{print $1 * 0.1}'`
      elif test ${gpair} = 'pr@rsds'; then
#        shdmax=`echo ${shdmax} | awk '{print $1 / 20.}'`
        cntmax=`echo ${cntmax} | awk '{print $1 / 3.}'`
      elif test ${gpair} = 'clt@cll' || test ${gpair} = 'clh@clm'; then
        cntcolor='red'
        cntfmt='%.1g'
      elif test ${gpair} = 'clh@clm'; then
        cntmax=`echo ${cntmax} | awk '{print $1}'`
      elif test ${gpair} = 'prw@huss'; then
        shdmax=`echo ${shdmax} | awk '{print $1*4}'`
      elif test ${gpair} = 'prls@prc'; then
        shdmax=`echo ${shdmax} | awk '{print $1/1.}'`
        cntmax=`echo ${cntmax} | awk '{print $1 * 0.5}'`
      elif test ${gpair} = 'hfls@hfss'; then
        cntmin='-50.'
        cntmax='50.'
      fi
  
      if test ${vark} = 'mean'; then
        graphvals=${shdstdn}','${cntstdn}':z|-1,y|-1:z|-1,y|-1:lat:pressure:'
        graphvals=${graphvals}${shdcbar}':fixsigc,'${cntcolor}':'${cntfmt}':'${shdmin}','
        graphvals=${graphvals}${shdmax}':'${cntmin}','${cntmax}',9:LMDZ+WRF|'${exp}
        graphvals=${graphvals}'|meridional|monthly|average|of|'${shdstdn}'|&|'
        graphvals=${graphvals}${cntstdn}':pdf:flip@y:None'
      else
        graphvals=${shdstdn}','${cntstdn}':y|-1,time|-1:y|-1,time|-1:time:XLAT:'
        graphvals=${graphvals}${shdcbar}':fixsigc,'${cntcolor}':'${cntfmt}':'${shdmin}','
        graphvals=${graphvals}${shdmax}':'${cntmin}','${cntmax}',9:WRF+LMDZ|'${exp}
        graphvals=${graphvals}'|mean|meridional|evolution|of|'${shdstdn}'|&|'
        graphvals=${graphvals}${cntstdn}':pdf:None:time|hours!since!1949-12-01|'
        graphvals=${graphvals}'exct,5,d|%d|date!([DD]):None'
      fi

      echo ${files}
      echo ${graphvals}
      echo ${shdstdn}'mean,'${cntstdn}'mean'

      if test ${vark} = 'sfc' || test ${vark} = 'diag' && 
       ! $(isin_list ${coup2D} ${gpair}); then
        python ${pyHOME}/drawing.py -f ${files} -o draw_2D_shad_cont_time -S ${graphvals} -v   \
          ${shdstdn}'mean,'${cntstdn}'mean'
        pyexc=$?
      else
        python ${pyHOME}/drawing.py -f ${files} -o draw_2D_shad_cont -S ${graphvals} -v        \
          ${shdstdn}'mean,'${cntstdn}'mean'
        pyexc=$?
      fi
      if test ${pyexc} -ne 0; then
        echo ${errormsg}
        echo "  "${main}": drawing.py fails!"
        exit
      else
        mv 2Dfields_shadow-contour.pdf ${ofold}/${exp}/${shdvar}mean_${cntvar}mean.pdf
        evince ${ofold}/${exp}/${shdvar}mean_${cntvar}mean.pdf &
      fi

      if $(isin_list ${coup2D} ${gpair}); then
        graphvals=${shdstdn}','${cntstdn}':south_north|-1,west_east|-1:'
        graphvals=${graphvals}'south_north|-1,west_east|-1:longitude:latitude:'
        graphvals=${graphvals}${shdcbar}':fixsigc,'${cntcolor}':'${cntfmt}':'${shdmin}','
        graphvals=${graphvals}${shdmax}':'${cntmin}','${cntmax}',9:LMDZ+WRF|'${exp}
        graphvals=${graphvals}'|monthly|average|of|'${shdstdn}'|&|'
        graphvals=${graphvals}${cntstdn}':pdf:flip@y:None'
        echo '  '${shdstdn}'mean,'${cntstdn}'mean'

        python ${pyHOME}/drawing.py -f ${files} -o draw_2D_shad_cont -S ${graphvals} -v        \
          ${shdstdn}'mean,'${cntstdn}'mean'
        pyexc=$?
        if test ${pyexc} -ne 0; then
          echo ${errormsg}
          echo "  "${main}": drawing.py fails!"
          exit
        else
          mv 2Dfields_shadow-contour.pdf ${ofold}/${exp}/${shdvar}tmean_${cntvar}tmean.pdf
          evince ${ofold}/${exp}/${shdvar}tmean_${cntvar}tmean.pdf &
        fi
      fi
#      exit

# end of couples
    done
# end of kinds
  done

done
fi

# 2D single plots
##
if ${single2Dplot}; then
  echo "2D single plots"
  for exp in ${exps}; do
    echo "  exp:"${exp}
    for tmean in ${ofold}/${exp}/*tmean.nc ;do
      fname=`basename ${tmean}`
      var=`ncdump -h ${tmean} | grep float | grep mean | tr '(' ' ' | awk '{print $2}'`
      if test ! ${var} = 'cltmean'; then
        shdvals=`python ${pyHOME}/drawing.py -o variable_values -S ${var} |          \
          grep all_values | awk '{print $3}'`
        cbar=`echo ${shdvals} | tr ',' ' ' | awk '{print $7}'`
        min=`echo ${shdvals} | tr ',' ' ' | awk '{print $3}'`
        max=`echo ${shdvals} | tr ',' ' ' | awk '{print $4}'`

        if test ${var} = 'prlsmean'; then xtrms='0,0.00015';
        elif test ${var} = 'prwmean'; then xtrms='0,40';
        elif test ${var} = 'psmean'; then xtrms='99000,102500';
        elif test ${var} = 'r2mean'; then xtrms='0,0.025';
        elif test ${var} = 'tasmean'; then xtrms='275,310';
        elif test ${var} = 'tmlamean'; then xtrms='260,310';
        elif test ${var} = 'uasmean'; then xtrms='-20,20';
        elif test ${var} = 'vasmean'; then xtrms='-20,20';
        else xtrms=${min}','${max}; fi

#        vals=${var}':XLONG|-1,XLAT|-1:XLONG:XLAT:'${cbar}':'${xtrms}':monthly|mean:pdf:'
#        vals=${vals}'None:None:true'
#        dims='south_north@XLAT,west_east@XLONG'
        python ${pyHOME}/nc_var.py -o WRF_CFlonlat_creation -S longitude,latitude    \
          -f ${tmean} -v XLONG,XLAT
        vals=${var}':longitude|-1,latitude|-1:longitude:latitude:'${cbar}':'${xtrms}':monthly|mean:pdf:'
        vals=${vals}'None:None:true'
        dims='south_north@latitude,west_east@longitude'
        python ${pyHOME}/drawing.py -o draw_2D_shad -S ${vals} -f ${tmean} -v ${var}
        pyexec=$?
        if test ${pyexec} -ne 0; then
          echo ${errormsg}
          echo "  "${main}": python fails!"
          exit
        else
          mv 2Dfields_shadow.pdf ${ofold}/${exp}/${var}tmean.pdf
#          evince ${ofold}/${exp}/${var}tmean.pdf &
#          exit
        fi
      fi
    done
  done
fi

# lon 2 vertical section
##
echo "lon 2 vertical section ..."
vars=`echo ${lon2DZsecvars} | tr ':' ' '`
pts=`echo ${lon2DZsecpts} | tr ':' ' '`

if ${lonZsec}; then
  for exp in ${exps}; do
    file=${ofold}/${exp}/vertical_interpolation_WRFp.nc
#    python ${pyHOME}/nc_var.py -o WRF_CFlonlat_creation -S longitude,latitude        \
#     -f ${file} -v lon,lat
## ALREADY done!
#    python ${pyHOME}/nc_var.py -o valmod -S lowthres@oper,0.,sumc,360. -f ${file}    \
#       -v lon
    if test ${exp} = 'AR40'; then labexp='wlmdza'
    else labexp='wlmdzb'; fi

    for pt in ${pts}; do
      echo "  pt: "${pt}
      vals='x:15|y:'${pt}'|time:23'
      lonval=`python ${pyHOME}/nc_var.py -o varout -f ${file} -S ${vals} -v lon |    \
        awk '{print $2}'`
      latval=`python ${pyHOME}/nc_var.py -o varout -f ${file} -S ${vals} -v lat |    \
        awk '{print $2}'`
      tval=`python ${pyHOME}/nc_var.py -o varout -f ${file} -S ${vals} -v time |     \
        awk '{print $2}'`
      for var in ${vars}; do
        echo "    var: "${var}

        shdvals=`python ${pyHOME}/drawing.py -o variable_values -S ${var} |          \
          grep all_values | awk '{print $3}'`
        cbar=`echo ${shdvals} | tr ',' ' ' | awk '{print $7}'`
        min=`echo ${shdvals} | tr ',' ' ' | awk '{print $3}'`
        max=`echo ${shdvals} | tr ',' ' ' | awk '{print $4}'`

# WRFt,U,V,WRFrh,WRFght
        if test ${var} = 'WRFght'; then xtrms='0,40000';
        elif test ${var} = 'WRFt'; then xtrms='200,300';
        elif test ${var} = 'U'; then xtrms='-40,40';
        elif test ${var} = 'V'; then xtrms='-20,20';
        else xtrms=${min}','${max}; fi

# Plotting
        vals=${var}':x|-1,y|'${pt}',z|-1,time|24:lon:pressure:'${cbar}':'
        vals=${vals}${xtrms}':'${labexp}'|vertical|longitudinal|section|('${latval}
        vals=${vals}'$^{\circ}$):pdf:flip@y:None:true'
        python ${pyHOME}/drawing.py -o draw_2D_shad -S ${vals} -f ${file} -v ${var}
        pyexec=$?
        if test ${pyexec} -ne 0; then
          echo ${errormsg}
          echo "  "${main}": python fails!"
          exit
        else
          mv 2Dfields_shadow.pdf ${ofold}/${exp}/${var}_lonZsec_${pt}pt.pdf
          evince ${ofold}/${exp}/${var}_lonZsec_${pt}pt.pdf &
        fi

        cntcolor='black'
        cntfmt='%g'

        if test ${var} = 'WRFght'; then cxtrms='0,40000';
        elif test ${var} = 'WRFt'; then cxtrms='200,300';
        elif test ${var} = 'U'; then cxtrms='-40,40';
        elif test ${var} = 'V'; then cxtrms='-20,20';
        else cxtrms=${min}','${max}; fi

        graphvals=${var}','${var}':x|-1,y|'${pt}',z|-1,time|24:'
        graphvals=${graphvals}'x|-1,y|'${pt}',z|-1,time|24:lon:pressure:'
        graphvals=${graphvals}${cbar}':fixsigc,'${cntcolor}':'${cntfmt}':'${xtrms}':'
        graphvals=${graphvals}${cxtrms}',15:'${labexp}
        graphvals=${graphvals}'|vertical|zonal|section|('${latval}'$^{\circ}$):pdf:'
        graphvals=${graphvals}'flip@y:None'

        python ${pyHOME}/drawing.py -o draw_2D_shad_cont -S ${graphvals}             \
          -f ${file},${file} -v ${var},${var}
        pyexec=$?
        if test ${pyexec} -ne 0; then
          echo ${errormsg}
          echo "  "${main}": python fails!"
          exit
        else
          mv 2Dfields_shadow-contour.pdf ${ofold}/${exp}/${var}_lonZsec-cnt_${pt}pt.pdf
          evince ${ofold}/${exp}/${var}_lonZsec-cnt_${pt}pt.pdf &
        fi

#        exit
      done
    done
  done
fi

echo "Computing differences"
exp1=`echo ${experiments} | tr ':' ' ' | awk '{print $1}'`
exp2=`echo ${experiments} | tr ':' ' ' | awk '{print $2}'`

diffks=`echo ${diffkinds} | tr ':' ' '`

if test ${differences}; then
  for kind in ${diffks}; do
    echo ${kind}
    case ${kind} in
      'diffZ')
        vars=`echo ${diffZvars} | tr ':' ' '`
        fileorigd='mean.nc'
      ;;
      'diffH')
        vars=`echo ${diffHvars} | tr ':' ' '`
        fileorigd='tmean.nc'
      ;;
    esac

    echo "  "${var}
    for var in ${vars}; do
      if test ! -f ${ofold}/${var}_diff.nc; then
#        cdo sub ${ofold}/${exp1}/${var}${fileorigd} ${ofold}/${exp2}/${var}${fileorigd} ${ofold}/${var}_diff.nc
        if test ${kind} = 'diffZ'; then
          values='pressure|lat@add|'${ofold}'/'${exp1}'/'${var}${fileorigd}'|'${var}'mean,'
          values=${values}'sub|'${ofold}'/'${exp2}'/'${var}${fileorigd}'|'${var}'mean'
          python ${pyHOME}/nc_var.py -o compute_opersvarsfiles -S ${values} -v ${var}
          pyexec=$?
          if test ${pyexec} -ne 0; then
            echo ${errormsg}
            echo "  "${main}": python fails!"
            exit
          fi
          mv opersvarsfiles_${var}.nc ${ofold}/${var}_diff.nc
        else
          values='longitude|latitude@add|'${ofold}'/'${exp1}'/'${var}${fileorigd}'|'${var}'mean,'
          values=${values}'sub|'${ofold}'/'${exp2}'/'${var}${fileorigd}'|'${var}'mean'
          python ${pyHOME}/nc_var.py -o compute_opersvarsfiles -S ${values} -v ${var}
          pyexec=$?
          if test ${pyexec} -ne 0; then
            echo ${errormsg}
            echo "  "${main}": python fails!"
            exit
          fi
          mv opersvarsfiles_${var}.nc ${ofold}/${var}_diff.nc
        fi
      fi
    done

    if test ${kind} = 'diffZ'; then
      coups=`echo ${graph_couples_mean} | tr ':' ' '`
      for coup in ${coups}; do
        shdvar=`echo ${coup} | tr '@' ' ' | awk '{print $1}'`
        cntvar=`echo ${coup} | tr '@' ' ' | awk '{print $2}'`

        echo "  couple: "${shdvar}'-'${cntvar}

        shdvals=`python ${pyHOME}/drawing.py -o variable_values -S ${shdvar} | grep all_values | awk '{print $3}'`
        cntvals=`python ${pyHOME}/drawing.py -o variable_values -S ${cntvar} | grep all_values | awk '{print $3}'`
        files=${ofold}'/'${shdvar}'_diff.nc,'${ofold}'/'${cntvar}'_diff.nc'

        shdstdn=`echo ${shdvals} | tr ',' ' ' | awk '{print $1}'`
        shdcbar=`echo ${shdvals} | tr ',' ' ' | awk '{print $7}'`
        shdmin=`echo ${shdvals} | tr ',' ' ' | awk '{print $3}'`
        shdmax=`echo ${shdvals} | tr ',' ' ' | awk '{print $4}'`

        cntstdn=`echo ${cntvals} | tr ',' ' ' | awk '{print $1}'`
        cntmin=`echo ${cntvals} | tr ',' ' ' | awk '{print $3}'`
        cntmax=`echo ${cntvals} | tr ',' ' ' | awk '{print $4}'`

        shdcbar='seismic'
        if test ${coup} = 'va@ua'; then
          shdmin=-4.
          shdmax=4.
          cntmin=-10.
          cntmax=10.
        elif test ${coup} = 'hus@ta'; then 
          shdmin=-0.2
          shdmax=0.2
          cntmin=-2.
          cntmax=2.
        fi

        cntcolor='black'
        cntfmt='%g'

        graphvals=${shdstdn}','${cntstdn}':x|-1,y|-1,x_2|-1:x|-1,y|-1,x_2|-1:lat:'
        graphvals=${graphvals}'pressure:'${shdcbar}':fixsigc,'${cntcolor}':'${cntfmt}':'
        graphvals=${graphvals}${shdmin}','${shdmax}':'${cntmin}','${cntmax}',9:'
        graphvals=${graphvals}${exp1}'-'${exp2}'|meridional|monthly|average|differences|of|'
        graphvals=${graphvals}${shdstdn}'|&|'${cntstdn}':pdf:flip@y:None'

        python ${pyHOME}/drawing.py -o draw_2D_shad_cont -S ${graphvals}             \
          -f ${files} -v ${shdvar},${cntvar}
        pyexec=$?
        if test ${pyexec} -ne 0; then
          echo ${errormsg}
          echo "  "${main}": python fails!"
          exit
        else
          mv 2Dfields_shadow-contour.pdf ${ofold}/${shdvar}-${cntvar}_diff.pdf
          evince ${ofold}/${shdvar}-${cntvar}_diff.pdf &
        fi
      done
#      exit
    else
      for var in ${vars}; do
        echo "  "${var}
        vals=`python ${pyHOME}/drawing.py -o variable_values -S ${var} | grep all_values | awk '{print $3}'`
        file=${ofold}/${var}_diff.nc

        stdn=`echo ${vals} | tr ',' ' ' | awk '{print $1}'`
        cbar='seismic'

        if test ${var} = 'uas'; then xtrms='-10.,10.';
        elif test ${var} = 'vas'; then xtrms='-5,5';
        elif test ${var} = 'ps'; then xtrms='-500,500';
        elif test ${var} = 'pr'; then xtrms='-0.001,0.001';
        else xtrms=${min}','${max}; fi

        vals=${var}':x|-1,y|-1:longitude:latitude:'${cbar}':'
        vals=${vals}${xtrms}':'${exp1}'-'${exp2}'|meridional|monthly|average|differences:'
        vals=${vals}'pdf:None:None:true'
        python ${pyHOME}/drawing.py -o draw_2D_shad -S ${vals} -f ${file} -v ${var}
        pyexec=$?
        if test ${pyexec} -ne 0; then
          echo ${errormsg}
          echo "  "${main}": python fails!"
          exit
        else
          mv 2Dfields_shadow.pdf ${ofold}/${var}_diff.pdf
          evince ${ofold}/${var}_diff.pdf &
        fi
#        exit
      done
    fi
  done
fi