#!/bin/ksh

# ICON
#
# ------------------------------------------
# Copyright (C) 2004-2026, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
# Contact information: icon-model.org
# See AUTHORS.TXT for a list of authors
# See LICENSES/ for license information
# SPDX-License-Identifier: BSD-3-Clause
# ------------------------------------------

#
# Cycle over periods (e.g. months, seasons, years),
# get monthly means and call multiple plot scripts:
#
# - difference maps against a control (compvar_ens.met)
# - difference maps against observations (climate_obs.met, Adrian Tompkins)
# - zonal means against control and ERA40 (zondia_seas_icon_batch.met,
#   Thomas Jung)
# - winds against ERA40 (wind_maps_clim.met, Peter Bechtold)
#
# AUTHORS:  Martin Koehler  April 2003 & April 2004
# UPDATED: Adrian Tompkins CY29R1 - new datasets
#                                   maps against ops
#                                   bug corrects
#                          CY29R2 - improved plots
#                                   new datasets
#                                   maps against seas
#          Martin Koehler  CY29R3 - parallel on linux cluster
#          Niles Wedi               & offline under SMS
#          Peter Bechtold  CY33R1 - verification against
#          Martin Koehler           ERA Interim
#          Peter Bechtold  CY38R1 - add wavenumber frequency statistics
#                                   adjust scripts for equal start times

#-------------------------------------------------------------------------#

#. $HOME/.kshrc               # for troubleshooting ecp

set -ex                      # exit if error, verbose

#use pgf90 # for compiling on desktop in offline version

export LM_LICENSE_FILE="/usr/local/apps/pgi/license.dat:$LM_LICENSE_FILE"


#############################################################
# Setup
#############################################################

#res=2.5
res=2.0

model='icon'
#model='ifs'

integer iperiod

# --- argument list processing

while getopts a:b:c:l:f:r:n:d:p:t:x:m:i:v:s:o:y:z:e:w:g:q: option
do
   case $option in
    a) EXPVER=$OPTARG;;      # experiment ID  of EXP (e.g. e8s7)
    b) EXPVER2=$OPTARG;;     # experiment ID2 of CONTROL (e.g. e9va or ERA40)
    c) CLASS=$OPTARG;;       # class
    d) DATESTR=$OPTARG;;     # multiple initial times (e.g. 19870501/19870502/19870503)
    f) FCLENGTH=$OPTARG;;    # forecast length [h]
    i) INITIME=$OPTARG;;     # initime
    l) NLEV=$OPTARG;;        # number of levels (e.g. 60 or 91)
    m) LCLIM=$OPTARG;;       # 1: plot against climatology
    s) LSIGT=$OPTARG;;       # 1: none 2: with statistical significance test
    o) LCOL=$OPTARG;;        # 1: colour plot, 2: b/w plot
    y) LARTIC=$OPTARG;;      # 1: standard     2: changes titles to a b c
    z) LDCYCLE=$OPTARG;;     # 1: plot diurnal cycle
    e) LMLDIF=$OPTARG;;      # 1: plot ml differences between 2 experiments
    q) LZONNOML=$OPTARG;;    # 1: do not make ml zonal plots (time consumimg)
    w) LWAVFREQ=$OPTARG;;    # 1: plot wave frequency diagrams (year only)
    g) LVDSPECT=$OPTARG;;    # 1: plot global vort div spectra  (year only)
    n) NUMBER=$OPTARG;;      # number of SEAS ensemble members
    p) PLOT_OPT=$OPTARG;;    # abc (e.g. 010) a=single month b=seasons c=whole exp mean
    r) FCRANGE=$OPTARG;;     # forecast range [d]
    t) TASK=$OPTARG;;        # task: getmeans, plot_dcycle, plot_2d, plot_wind, plot_z500, plot_obs_1-8, plot_clean
    v) CLASSAN=$OPTARG;;     # verifying Analysis ei or e4
    x) iperiod=$OPTARG;;     # period
   \?) errflg=1;;
   esac
done

if [[ $CLASSAN = "ei" ]] ; then
  eratext=ERAI
  expan=0001  # was 1112
else
  eratext=ERA40
  expan=1
fi

if [[ $errflg -ne 0 ]] ; then
  echo "incorrect usage!" >&2
  exit 1
fi

iperiod=iperiod-1            # ksh first index is 0

if [[ $EXPVER = $EXPVER2 ]] ; then
  DOPLOTDIFF=0
else
  DOPLOTDIFF=1               # do difference plots to 2nd exp. (1/0)
fi

# --- metview troubleshooting

#metview=metview_old
#metview=metview_dev
#metview=/usr/local/apps/Metview/metview4_new
metview=/usr/local/apps/Metview/4.7.0/bin/metview


case $TASK in
  "plot_2d")    # plot_2d  works with 4.6.0 nto 4.5.0 and not 4.7.0
    metview=/usr/local/apps/Metview/4.6.0/bin/metview
     ;;
  "plot_zon")   # plot_zon works with 4.5.0 not 4.6.0 and not 4.7.0
    metview=/usr/local/apps/Metview/4.5.0/bin/metview
     ;;
esac

#deia_grib_api=/home/ms/de/deia/grib_api # inlcudes fix add back in the MARS definitions for level reading
#                                        # this will be fixed in 1.12.3 by Gramalla
#dwd_grib_api=/home/ms/de/dwd/grib_api
#if [[ $TASK = "plot_2d" || $TASK = "plot_zon" || $TASK = "plot_obs_1" || $TASK = "plot_obs_3" || $TASK = "plot_obs_4" || $TASK = "plot_obs_5" ||  $TASK = "plot_obs_7" ||  $TASK = "plot_obs_10" ]] ; then
#  export METVIEW_EXTRA_GRIB_DEFINITION_PATH=${dwd_grib_api}/definitions.edzw-1.13.0-4:/usr/local/apps/grib_api/1.13.0/share/grib_api/definitions
#  export METVIEW_EXTRA_GRIB_DEFINITION_PATH=/usr/local/apps/grib_api/1.16.0/share/grib_api/definitions
#  ls
#else
#  export METVIEW_EXTRA_GRIB_DEFINITION_PATH=/usr/local/apps/grib_api/1.15.0/share/grib_api/definitions
#fi

#pgi_path="/usr/local/apps/nag/22/LP64"
#NAGLIB="-libnag.a"
#export MACRO_EXTRA_LIBS="-L$pgi_path/lib"
#export MACRO_EXTRA_F90_FLAGS=-Mfixed
#export GRIB_GRIBEX_MODE_ON=1
#export MACRO_EXTRA_LIBS="$MACRO_EXTRA_LIBS $NAGLIB"
#use 32bit version of NAGLIB for metview3
#NAGPATH="/usr/local/apps/nag/21/ILP32/lib"
#export MACRO_EXTRA_F90_FLAGS="-L$NAGPATH -lnag_nag -rpath $NAGPATH"

export MACRO_EXTRA_LIBS="-L/usr/local/apps/nag/22/LP64/lib -lnag"

#for interpolation of precipitation fields within metview - small values unit [m/s] not [mm/day]
export PRECIPITATION_THRESHOLD=0

# needed for some data files that have non-standard GRIB headers
export GRIBEX_CHECK='OFF'    # (climate_obs.met)

# read observational data from 'disc' on Linux cluster or 'ecfs'; default=disc
freadobs=disc                # (climate_obs.met)
obsdir=/scratch/ms/de/dei2/climplot.obs/
#obsdir=/perm/rd/pae/CLIMDATA/

if [ ! -d $obsdir ]; then
  print harddisk not available obsdata retrieval will be from ecfs
  freadobs=ecfs
  obsdir=/rdx/pa/data/
fi
print Obs read from $freadobs $obsdir

steps=24
ppfrq=24
dayl=24
# small planet scaling
rplradi=1.0
rcorioi=1.0

if [[ $steps -lt $ppfrq ]]; then
      steps=$ppfrq
fi
if [[ $steps -lt 1 ]]; then
      steps=1
fi
if [[ $dayl -lt 1 ]]; then
      dayl=1
fi

#############################################################
# time (get time information for processing of each period)
#############################################################

set -A nperiods `$metview -b dates.met $FCLENGTH $DATESTR $PLOT_OPT $dayl | grep nmon     | sed 's/nmon//'     | sed 'y/[],/   /'`
set -A yy_start `$metview -b dates.met $FCLENGTH $DATESTR $PLOT_OPT $dayl | grep yy_start | sed 's/yy_start//' | sed 'y/[],/   /'`
set -A mm_start `$metview -b dates.met $FCLENGTH $DATESTR $PLOT_OPT $dayl | grep mm_start | sed 's/mm_start//' | sed 'y/[],/   /'`
set -A mm_len   `$metview -b dates.met $FCLENGTH $DATESTR $PLOT_OPT $dayl | grep mm_len   | sed 's/mm_len//'   | sed 'y/[],/   /'`

# --- reverse order of periods: all-season-month
#dmk set -A yy_start `echo ${yy_start[*]} | awk '{for (i=NF; i>=1; i--) print $i}'`
#    set -A mm_start `echo ${mm_start[*]} | awk '{for (i=NF; i>=1; i--) print $i}'`
#xxx set -A mm_len   `echo ${mm_len[*]}   | awk '{for (i=NF; i>=1; i--) print $i}'`

set -A inidates `echo ${DATESTR} | sed 's/\// /g'`
set -A timestr  `echo ${INITIME} | sed 's/\// /g'`
set -A timestr2 `echo ${INITIME} | sed 's/\//,/g'`

if (( ${mm_start[${iperiod}]} < 10 )) ; then
  iniyearmonth=${yy_start[${iperiod}]}"0"${mm_start[${iperiod}]}
else
  iniyearmonth=${yy_start[${iperiod}]}${mm_start[${iperiod}]}
fi
yearmonth=${iniyearmonth}

yearmonthera=${yearmonth}

integer nens
integer i
date=$inidates
time=$timestr
if [[ ${NUMBER} = -99 ]] ; then
  nens=${#inidates[*]}
  # need to adjust time string for ANA retrieval if times equal
  i=1
  INITIMAN=${timestr[0]}
  while (( ${i} < $nens )) ; do
   if [[ ${timestr[i]} != ${timestr[i-1]} ]]; then
    INITIMAN=${INITIMAN}/${timestr[i]}
   fi
   i=i+1
  done
else
  nens=$NUMBER
  i=0
  num=1
  while (( ${i} < $nens )) ; do
    inidates=$inidates" "$date
    timestr=$timestr" "$time
    i=i+1
    num=${num}" "$i
  done
  set -A inidates $inidates
  set -A timestr $timestr
  set -A num $num
fi

#############################################################
# get data: monthly means from ERA and both experiments
#############################################################

if [[ $TASK = "getmeans" ]] ; then

#  rm -f ${PLOTDATA}/mmean_ERA40_*_alldates.grb
#  rm -f ${PLOTDATA}/mmean_${EXPVER}_*_alldates.grb
#  rm -f ${PLOTDATA}/mmean_${EXPVER2}_*_alldates.grb

  rm -f ${PLOTDATA}/mmean_ERA40_*.grb
  rm -f ${PLOTDATA}/mmean_${EXPVER}_*.grb
  rm -f ${PLOTDATA}/mmean_${EXPVER2}_*.grb

#############################################################
# --- read experiment 1 & 2 data (multiple ensemble members)
#############################################################

  i=0
  while (( ${i} < $nens )) ; do
    mm_len1=${mm_len[${iperiod}]}
    mm_len2=${mm_len1}
    iniyearmonth2=$iniyearmonth
    if [[ $mm_len2 -gt 12 ]]; then
      mm_len2="[12,$mm_len2]*"
      yy2a=$(print $inidates[i] |cut -c1-4)
      yy2b=$(print $iniyearmonth2 |cut -c1-4)
      mm2=$(print $iniyearmonth2 |cut -c5,6)
      iniyearmonth2="[$yy2b,$yy2a]*$mm2"
    fi
    mm_len1=${mm_len2}
    if [[ ${NUMBER} = -99 ]] ; then
      if [[ ${FCRANGE} = -99 ]] ; then
        MNFILE=mmean_${EXPVER}_*_ini${inidates[i]}_mon${iniyearmonth}_n${mm_len1}.grb
        MNFILE2=mmean_${EXPVER2}_*_ini${inidates[i]}_mon${iniyearmonth2}_n${mm_len2}.grb
        ECSUBDIR=${inidates[i]}${timestr[i]}/
        #ECSUBDIR="ec:icon/experiments/%SUITE%/%EXPNUM%/seasonal"
      else
        MNFILE=mmean_${EXPVER}_${CLASS}_day${FCRANGE}_*_mon${iniyearmonth}_n${mm_len1}.grb
        MNFILE2=mmean_${EXPVER2}_${CLASS}_day${FCRANGE}_*_mon${iniyearmonth2}_n${mm_len2}.grb
        ECSUBDIR=${CLASS}'/'
      fi
    else
        MNFILE=mmean_${EXPVER}_*_ini${inidates[i]}_mem${num[i]}_mon${iniyearmonth}_n${mm_len1}.grb
        MNFILE2=mmean_${EXPVER2}_*_ini${inidates[i]}_mem${num[i]}_mon${iniyearmonth2}_n${mm_len2}]}.grb
        ECSUBDIR=${inidates[i]}/
    fi
   #${ECFS_SYS_PATH}/ecp.p -o ec:/RDX/prepIFS/${EXPVER}/${ECSUBDIR}${MNFILE} ${PLOTDATA}
    ${ECFS_SYS_PATH}/ecp.p -o ${ECDIR_SEAS}/${MNFILE} ${PLOTDATA}
    if [[ $DOPLOTDIFF = 1 ]] ; then                                 # *=sfc pl ml
      echo " -- "
   #  ${ECFS_SYS_PATH}/ecp.p -o ec:/RDX/prepIFS/${EXPVER2}/${ECSUBDIR}${MNFILE2} ${PLOTDATA}
      ${ECFS_SYS_PATH}/ecp.p -o ${ECDIR2_SEAS}/${MNFILE2} ${PLOTDATA}
    fi
    i=i+1
  done

#############################################################
# --- read ERA data
#
# NOTE: * ERA hard-wired to 60 levels... not used if 91 level model version used...
#       * No CL type available yet for ERA Interim take ERA40 for climatology files
#
#############################################################

  # --- loop over months
  if [[ ${FCRANGE} = -99 ]] ; then

    imon=0
    while (( ${imon} < ${mm_len[${iperiod}]} )) ; do

#############################################################
# --- read ERA 40 data (time=00/06/12/18)
#############################################################

      if [[ $LCLIM = 1 ]] ; then
        cat > mars.requestclim << EOF
         retrieve,
            stream=mnth,
            time=$INITIMAN,
            levtype=sfc,
            expver=1,
            param=2T/2D/SD,
            class=e4,
            date=${yearmonthera}01,
            type=cl,
            grid=${res}/${res},
            target="${PLOTDATA}/mmean_ERA40_sfc_climate_ini${yearmonthera}01_mon${yearmonthera}.grb"
         retrieve,
            levelist=1/2/3/5/7/10/20/30/50/70/100/150/200/250/300/400/500/600/700/775/850/925/1000,
           # levelist=1/2/3/5/7/10/20/30/50/70/100/150/200/250/300/400/500/600/700/775/800/850/900/925/950/1000,
            stream=mnth,
            time=$INITIMAN,
            levtype=pl,
            expver=1,
            param=Z/T/R/U/V/W/Q,
            class=e4,
            date=${yearmonthera}01,
            type=cl,
            grid=${res}/${res},
            target="${PLOTDATA}/mmean_ERA40_pl_climate_ini${yearmonthera}01_mon${yearmonthera}.grb"
EOF
        if [[ $NLEV = 60 ]] ; then
         if [[ $LMLDIF != 0 || $LZONNOML = 0 || $LDCYCLE != 0  ]]; then
          cat >> mars.requestclim << EOF
           retrieve,
            levelist=1/TO/60,
            stream=mnth,
            time=$INITIMAN,
            levtype=ml,
            expver=1,
            param=T/U/V/W/Q/CC/CLWC/CIWC,
            class=e4,
            date=${yearmonthera}01,
            type=cl,
            grid=${res}/${res},
            target="${PLOTDATA}/mmean_ERA40_ml_climate_ini${yearmonthera}01_mon${yearmonthera}.grb"
EOF
         fi
        fi
        mars mars.requestclim
      fi

      cat > mars.requestan << EOF
        retrieve,
            stream=mnth,
            time=$INITIMAN,
            levtype=sfc,
            expver=$expan,
            param=2T/2D/SD,
            class=$CLASSAN,
            date=${yearmonthera}01,
            type=an,
            grid=${res}/${res},
            target="${PLOTDATA}/mmean_ERA40_sfc_ini${yearmonthera}01_mon${yearmonthera}.grb"
        retrieve,
            levelist=1/2/3/5/7/10/20/30/50/70/100/150/200/250/300/400/500/600/700/775/850/925/1000,
            #levelist=1/2/3/5/7/10/20/30/50/70/100/150/200/250/300/400/500/600/700/775/800/850/900/925/950/1000,
            stream=mnth,
            time=$INITIMAN,
            levtype=pl,
            expver=$expan,
            param=Z/T/R/U/V/W/Q/CC/CLWC/CIWC,
            class=$CLASSAN,
            date=${yearmonthera}01,
            type=an,
            grid=${res}/${res},
            target="${PLOTDATA}/mmean_ERA40_pl_ini${yearmonthera}01_mon${yearmonthera}.grb"
EOF
      if [[ $NLEV = 60 ]] ; then
        if [[ $LMLDIF != 0 || $LZONNOML = 0 || $LDCYCLE != 0  ]]; then
          cat >> mars.requestan << EOF
          retrieve,
     	    levelist=1/TO/60,
	        stream=mnth,
	        time=$INITIMAN,
	        levtype=ml,
	        expver=$expan,
	        param=T/U/V/W/Q/CC/CLWC/CIWC,
	        class=$CLASSAN,
	        date=${yearmonthera}01,
	        type=an,
            grid=${res}/${res},
            target="${PLOTDATA}/mmean_ERA40_ml_ini${yearmonthera}01_mon${yearmonthera}.grb"
EOF
       fi
      fi
      mars mars.requestan

###########################################################
# --- combine ERA data (cat)
###########################################################

      cat   ${PLOTDATA}/mmean_ERA40_sfc_ini*mon${yearmonthera}.grb \
        >>  ${PLOTDATA}/mmean_ERA40_sfc_alldates.grb
      cat   ${PLOTDATA}/mmean_ERA40_pl_ini*mon${yearmonthera}.grb \
        >>  ${PLOTDATA}/mmean_ERA40_pl_alldates.grb
      if [[ $LCLIM = 1 ]] ; then
        cat   ${PLOTDATA}/mmean_ERA40_sfc_climate_ini*mon${yearmonthera}.grb \
          >>  ${PLOTDATA}/mmean_ERA40_sfc_climate_alldates.grb
        cat   ${PLOTDATA}/mmean_ERA40_pl_climate_ini*mon${yearmonthera}.grb \
          >>  ${PLOTDATA}/mmean_ERA40_pl_climate_alldates.grb
      fi
      if [[ $LMLDIF != 0 || $LZONNOML = 0 || $LDCYCLE != 0  ]]; then
        if [[ ${NUMBER} = -99 && $NLEV = 60 ]] ; then
          cat   ${PLOTDATA}/mmean_ERA40_ml_ini*mon${yearmonthera}.grb \
            >>  ${PLOTDATA}/mmean_ERA40_ml_alldates.grb
          if [[ $LCLIM = 1 ]] ; then
            cat   ${PLOTDATA}/mmean_ERA40_ml_climate_ini*mon${yearmonthera}.grb \
             >>  ${PLOTDATA}/mmean_ERA40_ml_climate_alldates.grb
          fi
        fi
      fi

      rm -f ${PLOTDATA}/mmean_ERA40_*_ini*mon${yearmonthera}.grb

      month=`echo $yearmonthera | cut -c 5-6`
      if [[ $month = 12 ]] ; then
        let yearmonthera=yearmonthera+100-11
      else
        let yearmonthera=yearmonthera+1
      fi

      let imon=imon+1

    done  # end of loop over months (getting data)

  fi # fcrange

###########################################################
# --- combine model data (cat)
###########################################################

# this works for both SEAS and OPER since ini*mon wild card picks up
# both kinds of ensemble lables...

  mm_len1=${mm_len[${iperiod}]}
  mm_len2=${mm_len1}
  iniyearmonth2=$iniyearmonth
  if [[ $mm_len2 -gt 12 ]]; then
    mm_len2="[12,$mm_len2]*"
    yy2a=$(print $inidates[i] |cut -c1-4)
    yy2b=$(print $iniyearmonth2 |cut -c1-4)
    mm2=$(print $iniyearmonth2 |cut -c5,6)
    iniyearmonth2="[$yy2b,$yy2a]*$mm2"
  fi
  mm_len1=${mm_len2}
  if [[ ${FCRANGE} = -99 ]] ; then
    MNFILE=mmean_${EXPVER}_sfc_ini*mon${iniyearmonth}_n${mm_len1}.grb
    MNFILE2=mmean_${EXPVER2}_sfc_ini*mon${iniyearmonth2}_n${mm_len2}.grb
  else
    MNFILE=mmean_${EXPVER}_${CLASS}_day${FCRANGE}_*_mon${iniyearmonth}_n${mm_len1}.grb
    MNFILE2=mmean_${EXPVER2}_${CLASS}_day${FCRANGE}_*_mon${iniyearmonth2}_n${mm_len2}.grb
  fi

  cat   ${PLOTDATA}/${MNFILE} >>  ${PLOTDATA}/mmean_${EXPVER}_sfc_alldates.grb
  if [[ ${FCRANGE} = -99 ]] ; then
    cat   ${PLOTDATA}/mmean_${EXPVER}_pl_ini*mon${iniyearmonth}_n${mm_len1}.grb  \
      >>  ${PLOTDATA}/mmean_${EXPVER}_pl_alldates.grb
    if [[ $LMLDIF != 0 || $LZONNOML = 0 || $LDCYCLE != 0  ]]; then
      if [[ ${NUMBER} = -99 ]] ; then
        cat   ${PLOTDATA}/mmean_${EXPVER}_ml_ini*mon${iniyearmonth}_n${mm_len1}.grb  \
          >>  ${PLOTDATA}/mmean_${EXPVER}_ml_alldates.grb
      fi
    fi
  fi
#  rm -f ${PLOTDATA}/mmean_${EXPVER}_sfc_ini*mon${iniyearmonth}_n${mm_len1}.grb  \
#        ${PLOTDATA}/mmean_${EXPVER}_pl_ini*mon${iniyearmonth}_n${mm_len1}.grb   \
#        ${PLOTDATA}/mmean_${EXPVER}_ml_ini*mon${iniyearmonth}_n${mm_len1}.grb

  if [[ $DOPLOTDIFF = 1 ]] ; then
    cat   ${PLOTDATA}/${MNFILE2} >>  ${PLOTDATA}/mmean_${EXPVER2}_sfc_alldates.grb
    if [[ ${FCRANGE} = -99 ]] ; then
      cat   ${PLOTDATA}/mmean_${EXPVER2}_pl_ini*mon${iniyearmonth2}_n${mm_len2}.grb  \
        >>  ${PLOTDATA}/mmean_${EXPVER2}_pl_alldates.grb
      if [[ $LMLDIF != 0 || $LZONNOML = 0 || $LDCYCLE != 0  ]]; then
        if [[ ${NUMBER} = -99 ]] ; then
          cat   ${PLOTDATA}/mmean_${EXPVER2}_ml_ini*mon${iniyearmonth2}_n${mm_len2}.grb  \
            >>  ${PLOTDATA}/mmean_${EXPVER2}_ml_alldates.grb
        fi
      fi
    fi
    rm -f ${PLOTDATA}/mmean_${EXPVER2}_sfc_ini*mon${iniyearmonth2}_n${mm_len2}.grb  \
          ${PLOTDATA}/mmean_${EXPVER2}_pl_ini*mon${iniyearmonth2}_n${mm_len2}.grb   \
          ${PLOTDATA}/mmean_${EXPVER2}_ml_ini*mon${iniyearmonth2}_n${mm_len2}.grb
  fi

# fix to set centre=98 for pressure level data
#  mv  ${PLOTDATA}/mmean_${EXPVER}_pl_alldates.grb  ${PLOTDATA}/mmean_${EXPVER}_pl_alldates_centre78.grb
#   grib_set -s centre=98 ${PLOTDATA}/mmean_${EXPVER}_pl_alldates_centre78.grb \
#                         ${PLOTDATA}/mmean_${EXPVER}_pl_alldates.grb

fi


#############################################################
# 2D comparison plots
#############################################################

if [[ $TASK = "plot_2d" ]] ; then

  if [[ ${FCRANGE} = -99 ]] ; then

    interval=24                  # output interval [h] (also modify compvar_ens.met!)
    #interval=12                 # output interval [h] (also modify compvar_ens.met!)

    EXPVER3=${EXPVER2}           # experiment ID3 - to compare for extra variables

# --- variables and types (instantaneous or accumulated)
#     (SHORT LIST, 48)
#   set -A VARS    TCC LCC MCC HCC TP  CP  LSP SSHF SLHF BLH SSR STR TSR TTR 2T  2D  SD 10SI 10U 10V SST  SKT TCW TCWV T   T   T   T   T    T   T   T  TSRC TTRC TCLW TCIW CAPE
#   set -A TYPE    ins ins ins ins ins ins ins ins  ins  ins ins ins ins ins ins ins ins ins ins ins ins  ins ins ins ins ins ins ins ins  ins ins ins ins  ins  ins  ins  ins  ins
#   set -A NEXP     2   2   2   2   2   2   2   2    2    2   2   2   2   2   2   2   2   2   2   2   2    2   2   2   2   2   2   2   2    2   2   2   2    2    2    2    2    2
#   set -A LEV      0   0   0   0   0   0   0   0    0    0   0   0   0   0   0   0   0   0   0   0   0    0   0   0  925 850 700 500 300  200 100 50   0    0    0    0    0    0
#   set -A LEVTYPE sfc sfc sfc sfc sfc sfc sfc sfc  sfc  sfc sfc sfc sfc sfc sfc sfc sfc sfc sfc sfc sfc  sfc sfc sfc pl  pl  pl  pl  pl   pl  pl  pl  sfc  sfc  sfc  sfc  sfc  sfc
#   set -A STAT    avg avg avg avg avg avg avg avg  avg  avg avg avg avg avg avg avg avg avg avg avg avg  avg avg avg avg avg avg avg avg  avg avg avg avg  avg  avg  avg  avg  avg

# missing: BLH, TSRC, TTRC, optional: z500, sp

    set -A VARS    TCC LCC MCC HCC TP  CP  LSP SSHF SLHF SSR STR TSR TTR 2T  2D  SD 10SI 10U 10V SST SKT TCW TCWV Z   T   T   T   T   T   T   T   T   T  TCLW TCIW CAPE
    set -A TYPE    ins ins ins ins ins ins ins ins  ins  ins ins ins ins ins ins ins ins ins ins ins ins ins ins ins ins ins ins ins ins ins ins ins ins ins  ins  ins
    set -A NEXP     2   2   2   2   2   2   2   2    2    2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2    2    2
    set -A LEV      0   0   0   0   0   0   0   0    0    0   0   0   0   0   0   0   0   0   0   0   0   0   0  500 925 850 700 500 400 300 200 100 50   0    0    0
    set -A LEVTYPE sfc sfc sfc sfc sfc sfc sfc sfc  sfc  sfc sfc sfc sfc sfc sfc sfc sfc sfc sfc sfc sfc sfc sfc pl  pl  pl  pl  pl  pl  pl  pl  pl  pl  sfc  sfc  sfc
    set -A STAT    avg avg avg avg avg avg avg avg  avg  avg avg avg avg avg avg avg avg avg avg avg avg avg avg avg avg avg avg avg avg avg avg avg avg avg  avg  avg

# --- loop over variables

    NNN=0
    while (( ${NNN} < ${#VARS[*]} )) ; do
      echo '~~~~~~~~~~~~~~~~~~~~plotting '${VARS[${NNN}]}' ~~~~~~~~~~~~~~~~~~~at time:'
      date

# --- run metview

      if (( ${NEXP[${NNN}]} == 2 )) ; then
        $metview -b compvar_ens.met ${EXPVER} ${EXPVER2} ${DATESTR} ${iniyearmonth} \
          ${mm_len[${iperiod}]} ${interval} ${VARS[${NNN}]} \
          ${LEV[${NNN}]} ${LEVTYPE[${NNN}]} ${TYPE[${NNN}]} ${STAT[${NNN}]} ${res} \
          $LSIGT $LCOL $LARTIC $model
      else
        $metview -b compvar_ens.met ${EXPVER} ${EXPVER3} ${DATESTR} ${iniyearmonth} \
          ${mm_len[${iperiod}]} ${interval} ${VARS[${NNN}]} \
          ${LEV[${NNN}]} ${LEVTYPE[${NNN}]} ${TYPE[${NNN}]} ${STAT[${NNN}]} ${res} \
          $LSIGT $LCOL $LARTIC $model
      fi

      let NNN=NNN+1
    done     # end of loop over variables

  fi # fcrange

fi


#############################################################
# Wind maps against ERA
#############################################################

if [[ $TASK = "plot_wind" ]] ; then
  if [[ ${FCRANGE} = -99 ]] ; then
    $metview -b wind_maps_clim.met $EXPVER ${DATESTR} ${iniyearmonth} ${mm_len[${iperiod}]} \
                                   ${res} ${eratext}
  fi
fi


#############################################################
# Z500 against ERA
#############################################################

if [[ $TASK = "plot_z500" ]] ; then
  if [[ ${FCRANGE} = -99 ]] ; then
    $metview -b Z500_bias_era_mm.met $EXPVER ${DATESTR} ${iniyearmonth} ${mm_len[${iperiod}]} \
                                     ${res} ${eratext}
  fi
fi

#############################################################
# Diurnal cycle plots
#############################################################
if [[ $LDCYCLE = 1 ]] ; then
  if [[ $TASK = "plot_dcycle" ]] ; then
    if [[ ${FCRANGE} = -99 ]] ; then
     # plot internally limited to seasons YEAR DJF JJA
        $metview -b save_mean_diurnal_flux.met $EXPVER $CLASS ${DATESTR} ${INITIME} \
                                            ${iniyearmonth} ${mm_len[${iperiod}]}
# at moment in order to limit file numbers only seasons YEAR JJA DJF are plotted
        $metview -slog -b plot_amp_phase_clim.met $EXPVER ${DATESTR} ${iniyearmonth} ${mm_len[${iperiod}]} \
                                            ${freadobs} ${obsdir}
        rm -f ${PLOTDATA}/*.fld
    fi
  fi
fi

#############################################################
# Wavenumber Frequency plots
#############################################################
if [[ $LWAVFREQ = 1 ]] ; then
  if [[ $TASK = "plot_wavfreq" ]] ; then
    if [[ ${FCRANGE} = -99 ]] ; then
     # only plotted for YEAR
     if [[ ${mm_len[${iperiod}]} -eq 12 ]] ; then

        $metview -b marsretr_wavfreq.met $EXPVER fc $CLASS ${DATESTR} 365 ${INITIME} $steps

# convert from Grib to netcdf
       cd $PLOTDATA
       export MARS2NETCDF_CONFIG=${obsdir}utilities/wavfreq.cfnames
       for fil in tp_*.fld2  olr_*.fld2; do
         a=$( basename $fil .fld2 )
         mars -N<<EOF
         mars2netcdf,
         source=${a}.fld2,
         target=$a.nc,
         referencedate=AUTOMATIC
EOF
        done
        rm -f ${PLOTDATA}/*.fld2

# transform datestring to idl read format
        idldatestr=$( print ${DATESTR}|sed 's/\// /g' )
        set -A idldats $idldatestr
        integer ndi NDI
        NDI=${#idldats[*]}
        ndi=0
        IDLDATESTR=""
        while [ $ndi -lt NDI-1 ];do
          IDLDATESTR=${IDLDATESTR}"'${idldats[$ndi]}',"
          ndi=ndi+1
        done
        IDLDATESTR=${IDLDATESTR}"'${idldats[$ndi]}'"

# copy all idl files
        cp $SMSINCLUDE/*.pro $SMSINCLUDE/coastline*.dat  .
        chmod +x MAINgen.pro

# passing exp and dates to idl main
        ./MAINgen.pro $EXPVER ${IDLDATESTR} $rplradi
         sed  s/24.D/${steps}.D/ MAIN.pro >new.pro
         mv new.pro MAIN.pro

        use idl

# compile and run wave spectra using idl package
       idl << eof1
          .compile waveeq.pro
          .compile mapaxislab.pro
          .compile MAIN.pro
           MAIN
eof1
       filnad=_nens${NDI}_${iniyearmonth}_nmon${mm_len[${iperiod}]}.ps
       psselect -p2 olr*_wavfreq_*_nens_nmon12.ps olr_${EXPVER}_WavFreq$filnad
       psselect -p2 tp*_wavfreq_*_nens_nmon12.ps tp_${EXPVER}_WavFreq$filnad
# if u want to keep the f,k statistics as .nc fields comment line below
       rm -f tp*.wksts.* olr*.wksts.* tp*.mnavg.* olr*.mnavg.ps *.pro coastl*.dat
       rm -f olr*_wavfreq_*_nens_nmon12.ps tp*_wavfreq_*_nens_nmon12.ps
     fi
    fi
  fi
fi

#############################################################
# Global vorticity divergence spectra plots
#############################################################
if [[ $LVDSPECT = 1 ]] ; then
  if [[ $TASK = "plot_vdspectra" ]] ; then
    if [[ ${FCRANGE} = -99 ]] ; then
     # only computed for YEAR
     if [[ ${mm_len[${iperiod}]} -eq 12 ]] ; then
     #compute spectra, pass radius and lev
       lev=250
       cp $SMSINCLUDE/tcpack.F90  $SMSINCLUDE/tcpackens.F90 .
       chmod +x tcpack.F90 tcpackens.F90
       ./tcpack.F90 $rplradi
       use grib_api
       pgf90 -o spect tcpack.f90 $GRIB_API_INCLUDE $GRIB_API_LIB
       cp $SMSINCLUDE/*.pro .
       chmod +x rot_div_spectrum_plot.pro

       endstep=$((365*dayl))
       for lev in 250 50; do
         i=0
         rm -f out* fort.*
         while [ $i -lt $nens ];do
           cat<< EOF >dir
             retrieve,date=${inidates[i]},class=$CLASS,expver="${EXPVER}",levtype=pl,type=fc,
             level=$lev,param=D/VO,step=${endstep},time=${timestr[i]},target="targx1"
EOF
           mars dir
           ./spect
           cat< outx2 >> fort.7
           i=i+1
         done
       # extract spectral truncation from Grib file
         resol=$(grib_ls -p J  targx1 | awk -F" " '$1 ~ /^[0-9]/ &&  NF == 1 {print $0}' | uniq)
       # create ensemble mean
         ./tcpackens.F90 $resol $nens
         pgf90 tcpackens.f90
         a.out
         mv fort.8 curv.${lev}.rot.div.dat

       # plotting with idl
         ./rot_div_spectrum_plot.pro $lev $resol $EXPVER
         use idl
         idl << eof1
          .compile MAIN.pro
           MAIN
eof1
         filnad=_nens${nens}_${iniyearmonth}_nmon${mm_len[${iperiod}]}.ps
         mv *.ps $PLOTDATA/kespectra_${EXPVER}_P${lev}${filnad}
       done
     fi
    fi
  fi
fi

#############################################################
# Zonal mean plots
#############################################################

if [[ $TASK = "plot_zon" ]] ; then

# --- loop over plot types

  if [[ ${FCRANGE} = -99 ]] ; then
    if [[ ${NUMBER} = -99 ]] ; then
      if [[ $DOPLOTDIFF = 1 ]] ; then
       if [[ $LMLDIF = 1 ]] ; then
        set -A LEVTYPE pl      pl      ml      pl       ml
        set -A EXP1    $EXPVER $EXPVER $EXPVER $EXPVER  $EXPVER   # experiment
        set -A EXP2    ERA40   $EXPVER $EXPVER $EXPVER2 $EXPVER2  # control
        set -A nplot   5
       else
         if [[ $LZONNOML = 0 ]] ; then
           set -A LEVTYPE pl      pl      ml      pl
           set -A EXP1    $EXPVER $EXPVER $EXPVER $EXPVER         # experiment
           set -A EXP2    ERA40   $EXPVER $EXPVER $EXPVER2        # control
           set -A nplot   4
         else
           set -A LEVTYPE pl      pl      pl
           set -A EXP1    $EXPVER $EXPVER $EXPVER                 # experiment
           set -A EXP2    ERA40   $EXPVER $EXPVER2                # control
           set -A nplot   3
         fi
       fi
      else
        if [[ $LZONNOML = 0 ]] ; then
          set -A LEVTYPE pl      pl      ml
          set -A EXP1    $EXPVER $EXPVER $EXPVER                  # experiment
          set -A EXP2    ERA40   $EXPVER $EXPVER                  # control
          set -A nplot   3
        else
          set -A LEVTYPE pl      pl
          set -A EXP1    $EXPVER $EXPVER # experiment
          set -A EXP2    ERA40   $EXPVER # control
          set -A nplot   2
        fi
      fi
    else
      if [[ $DOPLOTDIFF = 1 ]] ; then
        set -A LEVTYPE pl      pl       pl
        set -A EXP1    $EXPVER $EXPVER  $EXPVER    # experiment
        set -A EXP2    ERA40   $EXPVER  $EXPVER2   # control
        set -A nplot   3
      else
        set -A LEVTYPE pl      pl
        set -A EXP1    $EXPVER $EXPVER                     # experiment
        set -A EXP2    ERA40   $EXPVER                     # control
        set -A nplot   2
      fi
    fi

# --- run metview: zonal plots

    iplot=0
    while (( ${iplot} < ${nplot} )) ; do

#   if (( iperiod+1 == nperiods )) ; then  # plot only all-month-mean zonal means

      $metview -b zondia_seas_icon_batch.met ${EXP2[${iplot}]} ${EXP1[${iplot}]} \
        ${inidates[0]} ${nens} ${iniyearmonth} ${mm_len[${iperiod}]} \
        ${LEVTYPE[${iplot}]} ${NLEV} avg ${res} ${eratext} ${model}

     #$metview -b zondia_seas_icon_batch.met ${EXP2[${iplot}]} ${EXP1[${iplot}]} \
     #  ${inidates[0]} ${nens} ${iniyearmonth} ${mm_len[${iperiod}]} \
     #  ${LEVTYPE[${iplot}]} ${NLEV} std ${res} ${eratext} ${model}

      let iplot=iplot+1
    done  # end of loop over zonal plots
  fi # fcrange

fi


#############################################################
# Maps against observations - 10 plots
#############################################################
# index of fields, loop here to save metview memory
# 1: lprecip
# 2: lradn
# 3: ltcwv
# 4: ltclw
# 5: ltciw
# 6: lcc
# 7: lwinds
# 8: lsflux
# 9: lciwc
# 10: l2d2tsd

TASK_A=`echo $TASK | cut -c 1-8`   # plot_obs_1 >>> plot_obs
if [[ $TASK_A = "plot_obs" ]] ; then

    plot=`echo $TASK | cut -c 10-`   # plot_obs_1 >>> 1

#   Now run twice...
    plotml=1
    if [[ $LMLDIF = 0 && $LZONNOML = 1 && $LDCYCLE = 0  ]]; then
      plotml=0
    fi
    $metview -b climate_obs.met $EXPVER $nens ${iniyearmonth} \
      ${mm_len[${iperiod}]} ${FCRANGE} 1 $plot 0  ${res} ${eratext} ${freadobs} ${obsdir} \
      $LSIGT $LCOL $LARTIC $plotml 24 $model
    if [[ $LCLIM = 1 ]] ; then
     $metview -b climate_obs.met $EXPVER $nens ${iniyearmonth} \
       ${mm_len[${iperiod}]} ${FCRANGE} 1 $plot 1  ${res} ERA40 ${freadobs} ${obsdir} \
       $LSIGT $LCOL $LARTIC $plotml 24 $model
    fi
fi


#############################################################
# clean up
#############################################################

if [[ $TASK = "plot_clean" ]] ; then

  if [[ 1 = 1 ]] ; then
    rm -f ${PLOTDATA}/mmean_ERA40_pl_alldates.grb
    rm -f ${PLOTDATA}/mmean_ERA40_ml_alldates.grb
    rm -f ${PLOTDATA}/mmean_${EXPVER}_sfc_alldates.grb
    rm -f ${PLOTDATA}/mmean_${EXPVER}_pl_alldates.grb
    rm -f ${PLOTDATA}/mmean_${EXPVER}_ml_alldates.grb
    rm -f ${PLOTDATA}/mmean_${EXPVER2}_sfc_alldates.grb
    rm -f ${PLOTDATA}/mmean_${EXPVER2}_pl_alldates.grb
    rm -f ${PLOTDATA}/mmean_${EXPVER2}_ml_alldates.grb
  fi

fi

exit
