Index: trunk/tools/model_graphics.bash
===================================================================
--- trunk/tools/model_graphics.bash	(revision 835)
+++ trunk/tools/model_graphics.bash	(revision 835)
@@ -0,0 +1,939 @@
+#!/bin/bash -x
+main='WRF+LMDZ_graphics.bash'
+errormsg='ERROR -- error -- ERROR -- error'
+HOMEpy=/home/lfita/etudes/WRF_LMDZ/svn/LMDZ_WRF/tools
+
+function variables_values(){
+## Function to provide values to plot the different variables
+#    return:
+#      1: var name
+#      2: std name
+#      3: minimum 
+#      4: maximum
+#      5: long name ('|' for spaces)
+#      6: units
+#      7: color palette (following: 
+#        http://matplotlib.org/1.3.1/examples/color/colormaps_reference.html)
+#
+#    $ variables_values 'WRFght'
+#    zg geopotential_height 0. 80000. geopotential|height m2s-2 rainbow
+
+  fname='variables_values'
+
+  variablesfile=${HOME}/bats/variables.inf
+
+  if test ! -f ${variablesfile}; then
+    echo ${errormsg}
+    echo "  "${fname}": variables file '"${variablesfile}"' does not exist !!!"
+  fi
+
+  case $1 in
+    'a_tht'|'LA_THT') varname='a_tht' ;;
+    'bils'|'LBILS') varname='bils' ;;
+    'clt'|'CLT'|'cldt'|'Total cloudiness') varname='clt' ;;
+    'cll'|'cldl'|'Low-level cloudiness') varname='cll' ;;
+    'clm'|'cldm'|'Mid-level cloudiness') varname='clm' ;;
+    'clh'|'cldh'|'High-level cloudiness') varname='clh' ;;
+    'dtcon'|'LDTCON') varname='dtcon' ;;
+    'evspsbl'|'QFX'|'evap'|'LEVAP') varname='evspsbl' ;;
+    'hfls'|'LH'|'flat'|'LFLAT') varname='hfls' ;;
+    'hfss'|'sens'|'LSENS') varname='hfss' ;;
+    'hus'|'WRFrh'|'LMDZrh') varname='hus' ;;
+    'huss'|'WRFrhs'|'LMDZrhs') varname='huss' ;;
+    'lat'|'XLAT') varname='lat' ;;
+    'lon'|'XLONG') varname='lon' ;;
+    'pr'|'RAINTOT'|'precip'|'Precip Totale liq+sol') varname='pr' ;;
+    'pracc'|'ACRAINTOT') varname='pracc' ;;
+    'prc'|'WRFprc'|'LPLUC'|'pluc') varname='prc' ;;
+    'pres'|'presnivs'|'pressure') varname='pres' ;;
+    'prls'|'WRFprls'|'LPLUL'|'plul') varname='prls' ;;
+    'prw'|'WRFprw') varname='prw' ;;
+    'ps'|'psfc'|'PSFC'|'psol'|'Surface Pressure') varname='ps' ;;
+    'q_th') varname='q_th' ;;
+    'r'|'QVAPOR') varname='r' ;;
+    'rsds'|'SWdnSFC'|'SWdn at surface'|'SWDOWN') varname='rsds' ;;
+    'rsdsacc') varname='rsdsacc' ;;
+    's_pblh'|'LS_PBLH') varname='zmla' ;;
+    's_pblt'|'LS_PBLT') varname='tmla' ;;
+    's_therm'|'LS_THERM') varname='s_therm' ;;
+    'ta'|'WRFt'|'temp'|'Air temperature') varname='ta' ;;
+    'tmlap'|'LS_PBLT'|'s_pblt') varname='tmla' ;;
+    'tas'|'T2'|'t2m'|'Temperature 2m') varname='tas' ;;
+    'time'|'time_counter') varname='time' ;;
+    'ua'|'vitu'|'U'|'Zonal wind') varname='ua' ;;
+    'uas'|'u10m'|'U10'|'Vent zonal 10m') varname='uas' ;;
+    'va'|'vitv'|'V'|'Meridional wind') varname='va' ;;
+    'vas'|'v10m'|'V10'|'Vent meridien 10m') varname='vas' ;;
+    'wake_h'|'LWAKE_H') varname='wakeh' ;;
+    'wake_s'|'LWAKE_S') varname='wakes' ;;
+    'wss'|'SPDUV') varname='wss' ;;
+    'xtime'|'XTIME') varname='xtime' ;;
+    'zg'|'WRFght'|'Geopotential height') varname='zg' ;;
+    'zmax_th'|'LZMAX_TH') varname='zmax_th' ;;
+    'zmla'|'LS_PBLH'|'s_pblh') varname='zmla' ;;
+    *)
+      echo $errormsg
+      echo "  "${fname}": variable '"$1"' not defined !!!"
+      exit
+    ;;
+  esac
+
+  varvals=`cat ${variablesfile} | grep ' '${varname}' ' | tr '@' ' '`
+  varvalsT=`cat ${variablesfile} | grep ' '${varname}' ' | tr ' ' '#'`
+  Lvarvals=`expr length ${varvalsT}'0'`
+
+  if test ${Lvarvals} -le 1; then
+    echo ${errormsg}
+    echo "  "${fname}": variable name '"${varname}"' not defined in file '"          \
+      ${variablesfile}"' !!!"
+    exit
+  fi
+  
+  echo ${varvals}
+}
+
+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
+    #######
+
+ifold='/bdd/PCER/workspace/lfita/etudes/WRF_LMDZ/WaquaL/WRF_LMDZ'
+#ofold='/san0/lflmd/etudes/WRF_LMDZ/WaquaL/LMDZ'
+ofold='/bdd/PCER/workspace/lfita/etudes/WRF_LMDZ/WaquaL/WRF_LMDZ'
+
+experiments='AR40:NPv31'
+#varkinds='mean:sfc:diag'
+varkinds='sfc'
+
+mean_variables='dim2mean|WRFt:dim2mean|U:dim2mean|V:dim2mean|WRFght:dim2mean|WRFrh'
+
+# Variables to compute longitudinal means for a Hovmoeller
+sfc_variables='dim2mean|T2:dim2mean|U10:dim2mean|V10:dim2mean|PSFC:dim2mean|RAINTOT:'
+sfc_variables=${sfc_variables}'dim2mean|SWDOWN:dim2mean|WRFprw:dim2mean|WRFrhs:'
+sfc_variables=${sfc_variables}'dim2mean|HFX:dim2mean|LH:dim2mean|Q2:'
+sfc_variables=${sfc_variables}'dim2mean|PBLH:dim2mean|QFX:'
+sfc_variables='dim2mean|RAINTOT'
+sfc_variables_AR40=${sfc_variables}
+sfc_variables_NPv31=${sfc_variables}':dim2mean|LWAKE_S:dim2mean|LWAKE_H:dim2mean|LS_THERM:'
+sfc_variables_NPv31=${sfc_variables_NPv31}'dim2mean|LZMAX_TH'
+
+# Variables to get from the diagnostics
+diag_variables='dim2mean|clt:dim2mean|cllmh:'
+diag_variables=${diag_variables}'dim2mean|WRFprc:dim2mean|WRFprls:'
+#diag_variables=${diag_variables}'dim2mean|WRFbils'
+
+# Variables to compute temporal means on a 2D map
+#var2Dtmeans_AR40='WRFprc:WRFprls:PBLH:QFX:HFX:LH:WRFbils'
+var2Dtmeans_AR40='WRFprc:WRFprls:PBLH:QFX:HFX:LH'
+var2Dtmeans_AR40='WRFprc:WRFprls:PBLH:QFX:HFX:LH:WRFbils:RAINTOT:U10:V10:T2:'
+var2Dtmeans_AR40=${var2Dtmeans_AR40}'PSFC:clt:cll:clm:clh:cllmh:WRFprw:Q2:pr:prw'
+
+var2Dtmeans_NPv31=${var2Dtmeans_AR40}':LWAKE_S:LWAKE_H:LS_THERM:LZMAX_TH'
+
+#inst_variables='mclt_lmdz:pr_lmdz:mwss'
+graph_couples_mean='hus@ta:va@ua'
+graph_couples_sfc='tas@ps:uas@vas:pr@rsds:prw@huss:prls@prc:'
+#graph_couples_sfc=${graph_couples_sfc}'hfls@hfss:evspsbl@bils'
+graph_couples_sfc=${graph_couples_sfc}'hfls@hfss'
+graph_couples_sfc_AR40=${graph_couples_sfc}
+graph_couples_sfc_NPv31=${graph_couples_sfc}':zmaxth@stherm:wakeh@wakes:zmaxth@stherm'
+graph_couples_diag='prls@prc:clt@cll:clh@clm'
+#graph_couples_z='dtcon@pluc'
+
+#graph_couples_diag='prls@prc'
+
+lon2DZsecvars='WRFt:U:V:WRFrh:WRFght'
+lon2DZsecpts=9:15:21
+
+# Differences between experiments
+diffkinds='diffZ:diffH'
+diffZvars='ua:va:ta:hus'
+diffHvars='uas:vas:ps:pr'
+
+filenmean='vertical_interpolation_WRFp.nc'
+#ifoldmean='/media/data1/etudes/WRF_LMDZ/WaquaL/WRF_LMDZ'
+ifoldmean='/bdd/PCER/workspace/lfita/etudes/WRF_LMDZ/WaquaL/WRF_LMDZ'
+
+filensfc='wrfout_d01_1980-03-01_00:00:00'
+#ifoldsfc='/media/data1/etudes/WRF_LMDZ/WaquaL/WRF_LMDZ'
+ifoldsfc='/bdd/PCER/workspace/lfita/etudes/WRF_LMDZ/WaquaL/WRF_LMDZ'
+
+filendiag='diagnostics.nc'
+ifolddiag='/bdd/PCER/workspace/lfita/etudes/WRF_LMDZ/WaquaL/WRF_LMDZ'
+
+timeval='tstep'
+zval='null'
+dimval='lon,lat,pressure,time'
+kindval='pdf'
+
+compute=0
+computetlon=false
+plot=false
+single2Dplot=false
+lonZsec=false
+differences=true
+
+plevels='100000.,97500.,95000.,92500.,90000.,85000.,80000.,75000.,70000.,'
+plevels=${plevels}'65000.,60000.,55000.,50000.,45000.,40000.,35000.,30000.,'
+plevels=${plevels}'25000.,20000.,15000.,10000.,5000.,2500.,1000.,'
+plevels=${plevels}'500.,250.'
+
+# Couples to draw on a 2D map
+coup2D_AR40='prls@prc::zmla@tmla:hfls@hfss:evspsbl@bils'
+coup2D_NPv31=${coup2D_AR40}':wakeh@wakes:zmaxth@stherm'
+
+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
