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

main='model_graphics.bash'
errormsg='ERROR -- error -- ERROR -- error'
warnmsg='ERROR -- error -- ERROR -- error'

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 isin_list() {
# 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
}

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

# Uploading environment
uploadvars model_graphics.dat

export 
exit


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

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

combosfile=${HOMEpy}/diagnostics.inf

errmsg='ERROR -- error -- ERROR -- error'
####### ###### ##### #### ### ## #

exps=`echo ${experiments} | tr ':' ' '`
varks=`echo ${varkinds} | tr ':' ' '`

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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/nc_var.py -o WRF_CFxtime_creation -f ${ifoldsfc}/${exp}/${filensfc} \
            -S 19491201000000,hours -v time
          python ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/diagnostics.py -d ${values} -f ${ifoldsfc}/${exp}/${filensfc} -v \
            ${varcombos}
          if test $? -ne 0; then
            echo ${errormsg}
            echo "  python failed!!"
            echo python ${HOMEpy}/diagnostics.py -d ${values} -f ${ifoldsfc}/${exp}/${filensfc}\
              -v ${varcombos}
            exit
          fi
          python ${HOMEpy}/nc_var.py -o WRF_CFxtime_creation -f diagnostics.nc                 \
            -S 19491201000000,hours -v time
           python ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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
        echo "HereLluis"${var}
# 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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/drawing.py -o variable_values -S ${shdvar} | grep all_values | awk '{print $3}'`
      cntvals=`python ${HOMEpy}/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 ${HOMEpy}/drawing.py -f ${files} -o draw_2D_shad_cont_time -S ${graphvals} -v   \
          ${shdstdn}'mean,'${cntstdn}'mean'
        pyexc=$?
      else
        python ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/nc_var.py -o WRF_CFlonlat_creation -S longitude,latitude        \
#     -f ${file} -v lon,lat
## ALREADY done!
#    python ${HOMEpy}/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 ${HOMEpy}/nc_var.py -o varout -f ${file} -S ${vals} -v lon |    \
        awk '{print $2}'`
      latval=`python ${HOMEpy}/nc_var.py -o varout -f ${file} -S ${vals} -v lat |    \
        awk '{print $2}'`
      tval=`python ${HOMEpy}/nc_var.py -o varout -f ${file} -S ${vals} -v time |     \
        awk '{print $2}'`
      for var in ${vars}; do
        echo "    var: "${var}

        shdvals=`python ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/drawing.py -o variable_values -S ${shdvar} | grep all_values | awk '{print $3}'`
        cntvals=`python ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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 ${HOMEpy}/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
