###############################################################################
#    This file gathers a set of bash functions that rely on cdftools to       #
#                                                                             #
# reduce_mmo                                                                  #
# get_diagsMMO                                                                #
# get_nemovar                                                                 #
# get_glorys                                                                  #
# clean_diagsMMO                                                              #
# vertmeansal                                                                 #
# heat_sal_mxl                                                                #
# ohc_specified_layer                                                         #
# moc                                                                         #
# convection                                                                  #
# psi                                                                         #
# gyres                                                                       #
# area_moc                                                                    #
# max_moc                                                                     #
# siasiesiv                                                                   #
# ohc                                                                         #
# cutsection                                                                  #
# interp3d                                                                    #
# setminmax                                                                   #
# concat                                                                      #
# gather_memb                                                                 #
#                                                                             #
#      Those functions would never have seen the day without Hui Du,          #
#                     usually referred to as Super-Hui.                       #
#                                                                             #
#   He made a crucial work to develop what ended up below in the functions    #
#   that computes the sea ice extent, sea ice area, ocean heat content and    #
#   meridional overturning streamfunction.                                    #
#   Especially, he developped new options from the cdftools sources to be     #
#   able to compute the heat content in different basins.                     #
#                                                                             #
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#               You want to make available a new diagnostic ?                 @
#                                                                             @
#  1) Get an MMO tar files from any experiment on cfunas                      @
#  2) Write a bash function that works on a grid_T or grid_U or grid_V or     @
#     grid_W file or a combination from this MMO file                         @
#  --> You can test your function by defining the CON_FILES and NEMOVERSION   @
#      variables and by sourcing the current file, the meshmasks will be      @
#      available after sourcing, remember to source again after any           @
#      modification of your function                                          @
#  --> Your function should work on input files of any resolution             @
#      ORCA1/ORCA025/ORCA2                                                    @
#  --> Your function should work on input files of any time length            @
#  --> The output file should contain a proper time axis that you can copy    @
#      from your input file                                                   @
#  --> The output file should be at most a 3d field including the time        @
#      dimension                                                              @
#  3) Write a short description of your function, and add its name to the     @
#     list above                                                              @
#  4) Go the the ocean_pp.sh script to add a call to your function            @
#                                                                             @
#                      Any doubt ---> vguemas@ic3.cat                         @
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#                         Link constant file for co                           #
###############################################################################

ln -sf ${CON_FILES}/mesh_mask_nemo.${NEMOVERSION}.nc mesh_hgr.nc
ln -sf ${CON_FILES}/mesh_mask_nemo.${NEMOVERSION}.nc mesh_zgr.nc
ln -sf ${CON_FILES}/mesh_mask_nemo.${NEMOVERSION}.nc mask.nc
ln -sf ${CON_FILES}/new_maskglo.${NEMOVERSION}.nc new_maskglo.nc
if [[ ! -f mask.nc ]] ; then 
  echo "No configuration files for cdftools"
  exit 
fi 

###############################################################################
#       Reduced number of variables in diag files to save disk space          #
#                                                                             #
# $1 : input grid_T file name                                                 #
# $2 : input icemod file name                                                 #
# $3 : suffix of output files with nc extension                               #
#                                                                             #
#         Created in February 2012      Author : vguemas@ic3.cat              #
#         May 2014 : Compatibility with PA changes to oce output - Virginie   #
###############################################################################

function reduce_mmo {
ncks -O -v sosstsst,sosaline,somixhgt,somxl010,sossheig $1 oce_${3} 
typeset var lstvars=`cdo showvar $2`
if [[ ${lstvars/ileadfra} != ${lstvars} ]]  ; then
  ncks -O -v isnowthi,iicethic,ileadfra,iicetemp $2 ice_${3} 
else 
  ncks -O -v isnowthi,iicethic,iiceconc,iicetemp $2 ice_${3}
fi
ncks -O -v votemper $1 t3d_${3} 
}
###############################################################################
#                    Copy diags or MMO files from cfunas                      #
#                                                                             #
# $1 : starting date                                                          #
# $2 : expid                                                                  #
# $3 : member                                                                 #
# $4 : starting leadtime                                                      # 
# $5 : end leadtime                                                           #
# $6 : chunk length in month                                                  #
# $7 : nemo/ecearth                                                           #
# $8 : diags/MMO                                                              #
# $9 : storage frequency (daily/monthly)                                      #
# $10 : list of files extracted from the original tarballs
#                                                                             #
#         Created in May 2012           Author : vguemas@ic3.cat              #
#         Option 10: June 2013          isabel.andreu-burillo@ic3.cat
###############################################################################

function get_diagsMMO {
typeset var yyyy0=`echo $1|cut -c1-4`
typeset var mm0=`echo $1|cut -c5-6`
if [ -z "${10}" ] ; then
  typeset var lstypes="grid_T grid_U grid_V grid_W icemod"
else
  typeset var lstypes=${10}
fi

typeset var jt
typeset var year1
typeset var year2
typeset var mon1
typeset var mon2
for jt in $(seq $4 $6 $5) ; do
  year1=$(($yyyy0+(10#$mm0+$jt-2)/12)) 
  mon1=$(((10#$mm0+$jt-2)%12+1))
  year2=$(($yyyy0+(10#$mm0+$jt+$6-3)/12))
  mon2=$(((10#$mm0+$jt+$6-3)%12+1)) 
  cp /cfunas/exp/$7/$2/$1/fc$3/outputs/$8_$2_$1_fc$3_${year1}$(printf "%02d" $mon1)01-${year2}$(printf "%02d" $mon2)*.tar .
 if [[ "$8" == "MMO" ]]; then 
  for filetype in $lstypes; do
   tar --wildcards -xvf $8_$2_$1_fc$3_${year1}$(printf "%02d" $mon1)01-${year2}$(printf "%02d" $mon2)*.tar "*${freqkeep}*${filetype}*"
  done
 else
   tar --wildcards -xvf $8_$2_$1_fc$3_${year1}$(printf "%02d" $mon1)01-${year2}$(printf "%02d" $mon2)*.tar
 fi
  rm -f $8_$2_$1_fc$3_${year1}$(printf "%02d" $mon1)01-${year2}$(printf "%02d" $mon2)*.tar
  if [[ `ls *.gz` != '' ]] ; then gunzip -f *.gz ; fi
done

typeset var listroots
case $8 in 
  'diags' ) 
    listroots="t3d heat_sal_mxl ice moc psi sal_0-300m sal_300-5400m" 
    if [[ `ls sstsimld*`  != '' ]] ; then
      listroots=$listroots" sstsimld"
    elif [[ `ls sstsssmld*`  != '' ]] ; then
      listroots=$listroots" sstsssmld"
    else
      listroots=$listroots" oce"
    fi 
    ;;
  'MMO' ) listroots=$lstypes ;;
esac
case $9 in
  'daily') freqexcl1='1m' ; freqexcl2='MM' ;;
  'monthly' ) freqexcl1='1d' ; freqexcl2='DD' ;;
  *) freqexcl1='1d' ; freqexcl2='DD' ;;
esac

function concat_startdate {
  typeset var root
  typeset var lstfiles
  for root in ${listroots[@]} ; do
    if [[ `ls *${root}*` != '' ]] && [[ `ls *${root}*` != ${root}_1m.nc ]] ; then
     if [[ "$8" == "MMO" ]] ; then
      lstfiles=`ls *${root}* | grep -v "${root}_$2_$1_fc" | grep -v "${freqexcl1}" | grep -v "${freqexcl2}" | grep -v "km" `
     else
      lstfiles=`ls ${root}* | grep -v "${root}_$2_$1_fc" | grep -v "${freqexcl1}" | grep -v "${freqexcl2}" | grep -v "km" `
     fi

     file1=`echo "${lstfiles}" | tail -1`
     cdo_version=`cdo -V &> ff ; grep Climate ff | cut -d \  -f 5`
     rm ff
     #test on cdo version: if >1.5.6, remove valid_min/max attributes to avoid values out of that range to be replaced by NaN
     if [[ "$cdo_version" = "`echo -e "$cdo_version\n1.5.6" | sort -V | head -n1`" ]] ; then
      if [[ $root == 'grid_T' || $root == 't3d' ]] ; then
       for file in $lstfiles ; do
        ncatted -O -a valid_max,votemper,d,, $file $file
        ncatted -O -a valid_min,votemper,d,, $file $file
       done
      fi
      if [[ $root == 'heat_sal_mxl' ]] ; then
       for file in $lstfiles ; do
        ncatted -O -a valid_max,somxlheatc,d,, $file $file
        ncatted -O -a valid_min,somxlheatc,d,, $file $file
       done
      fi
     fi

     outfile=${root}_$2_$1_fc$3_$(($yyyy0+(10#$mm0+$4-2)/12))$(printf "%02d" $(((10#$mm0+$4-2)%12+1)))_${year2}$(printf "%02d" $mon2).nc
     typeset var lstvars=`cdo showvar $file1`
     if [[ ${lstvars/iicenflx} != ${lstvars} ]] ; then for file in $lstfiles ; do ncks -O -x -v iicenflx $file $file ; done ; fi
     cdo mergetime $lstfiles ${outfile}
     timevar=`ncdump -h ${outfile} | grep UNLIMITED | awk '{print $1}'`
     if [[ $timevar == 'time_counter' ]] ; then ncrename -v time_counter,time -d time_counter,time ${outfile} ; fi
     if [[ $root == 'moc' ]] ; then 
      lstdims=`ncdump -h ${outfile} | awk /dimensions:/,/variables:/ | grep -v dimensions: | grep -v variables: | awk '{print $1}'`
      if [[ ${lstdims/gsize} != ${lstdims} ]] ; then
       ncrename -d gsize,y ${outfile} 
      fi
      lenx=`ncdump -h ${outfile} | grep 'x =' | head -1 | awk '{print $3}'`
      if [[ $lenx > 1 ]] ; then
       if [[ ${lstvars/nav_lon} != ${lstvars} ]] ; then
         ncks -O -x -v nav_lon,nav_lat ${outfile} ${outfile}
       fi
       ncrename -d x,y ${outfile}
      fi
      ncks -A -v nav_lon,nav_lat `echo $lstfiles | awk '{print $1}' ` ${outfile}
     fi
     rm -f $lstfiles 
     if [[ $root == 'sstsimld' || $root == 'sstsssmld' ]] ; then mv ${outfile} oce_$2_$1_fc$3_$(($yyyy0+(10#$mm0+$4-2)/12))$(printf "%02d" $(((10#$mm0+$4-2)%12+1)))_${year2}$(printf "%02d" $mon2).nc ; fi
    fi
  done
}

concat_startdate $1 $2 $3 $4 $5 $6 $7 $8 

# These lines aim at concatenating the daily means as well and computing the monthly means from these daily means
if [[ $9 == 'monthly' ]] ; then
  freqexcl1='1m' ; freqexcl2='MM'
  for root in ${listroots[@]} ; do
    outfile=${root}_$2_$1_fc$3_$(($yyyy0+(10#$mm0+$4-2)/12))$(printf "%02d" $(((10#$mm0+$4-2)%12+1)))_${year2}$(printf "%02d" $mon2).nc
    if [[ -e $outfile ]] ; then
      mv $outfile ${root}_1m.nc
    fi
  done
  concat_startdate $1 $2 $3 $4 $5 $6 $7 $8
  for root in ${listroots[@]} ; do
    outfile=${root}_$2_$1_fc$3_$(($yyyy0+(10#$mm0+$4-2)/12))$(printf "%02d" $(((10#$mm0+$4-2)%12+1)))_${year2}$(printf "%02d" $mon2).nc
    if [[ -e $outfile ]] ; then
      cdo monmean $outfile ${root}_daily2monthly.nc
      rm -f $outfile
      if [[ -e ${root}_1m.nc ]] ; then
        mv ${root}_1m.nc $outfile
        ncks -A ${root}_daily2monthly.nc $outfile
        rm -f ${root}_daily2monthly.nc
      else
        mv ${root}_daily2monthly.nc $outfile
      fi
    else
      if [[ -e ${root}_1m.nc ]] ; then
        mv ${root}_1m.nc $outfile
      fi
    fi
  done
fi

rm -f *${freqexcl1}* *${freqexcl2}*

}
###############################################################################
#                         Copy NEMOVAR files from cfunas                      #
#                                                                             #
# $1 : expid                                                                  #
# $2 : member                                                                 #
# $3 : start year                                                             # 
# $4 : end year                                                               #
# $5 : start month                                                            #
# $6 : end month                                                              #
# $7 : list of files extracted from the original tarballs
#                                                                             #
#         Created in May 2012           Author : vguemas@ic3.cat              #
#         Modified: June 2013           isabel.andreu-burillo@ic3.cat         #
###############################################################################

function get_nemovar {

if [ -z "$5" ] ; then
  typeset var moni=9
else
  typeset var moni=$5
fi

if [ -z "$5" ] ; then
  typeset var monf=8
else
  typeset var monf=$6
fi

typeset var path
typeset var yearf
case $1 in
  'nemovar_s4') path=/cfunas/exp/ECMWF/NEMOVAR_S4/outputs/fc$2/s4 ;;
  'nemovar_combine') path=/cfunas/exp/ECMWF/NEMOVAR_COMBINE/outputs/opa0/fa9p_1m ;;
esac
typeset var year
typeset var mon
for year in $(seq $3 $4) ; do
  case $year in
    $3) mona=${moni} ;;
    *) mona=1 ;;
  esac
  case $year in
    $4) monb=${monf} ;;
    *) monb=12 ;;
  esac   
  for mon in $(seq $mona $monb); do
    cp ${path}_fc$2_${year}$(printf "%02d" $mon)*.gz .
  done
done
gunzip -f *.gz

typeset var listroots=${7}
typeset var root
typeset var lstfiles
typeset var ntimes
typeset var jt
for root in ${listroots[@]} ; do
  lstfiles=`ls *fc${2}*${root}* | grep -v "${root}_$1_195709_fc$2_${3}09_${4}$(printf "%02d" $monb).nc"`
  ncrcat -O -x -v vorbiasp $lstfiles tmp_${root}.nc 
  cdo settaxis,${3}-$(printf "%02d" $moni)-15,12:00,1mon tmp_${root}.nc ${root}_$1_19570901_fc$2_${3}$(printf "%02d" $moni)_${4}$(printf "%02d" $monf).nc
  rm -f $lstfiles tmp_${root}.nc
done
}
###############################################################################
#                         Copy GLORYS files from cfunas                       #
#                                                                             #
# $1 : start year                                                             # 
# $2 : end year                                                               #
# $3 : start month                                                            #
# $4 : end month                                                              #
#                                                                             #
#         Created in June 2013          Author : vguemas@ic3.cat              #
###############################################################################

function get_glorys {
typeset var path=/cfunas/exp/MERCATOR/GLORYS2V1/outputs/ORCA1L46   #ORCA025L75_glorys
typeset var lstfiles=""
for year in $(seq $1 $2) ; do
  cp ${path}/vosaline_${year}.nc .
  cp ${path}/votemper_${year}.nc .
  ncks -A vosaline_${year}.nc votemper_${year}.nc
  rm -f vosaline_${year}.nc
  lstfiles=${lstfiles}" "votemper_${year}.nc
done
cdo cat ${lstfiles} tmp.nc
cdo settaxis,${1}-01-15,12:00,1mon tmp.nc tmp2.nc
cdo seldate,${1}-$(printf "%02d" $3)-00,${2}-$(printf "%02d" $4)-31 tmp2.nc grid_T_glorys2v1_19930101_fc0_${1}$(printf "%02d" $3)_${2}$(printf "%02d" $4).nc
rm -f ${lstfiles} tmp.nc tmp2.nc
ncks -O -x -v nav_lon,nav_lat,x_2,y_2 grid_T_glorys2v1_19930101_fc0_${1}$(printf "%02d" $3)_${2}$(printf "%02d" $4).nc grid_T_glorys2v1_19930101_fc0_${1}$(printf "%02d" $3)_${2}$(printf "%02d" $4).nc
ncks -A -v nav_lon,nav_lat mesh_hgr.nc grid_T_glorys2v1_19930101_fc0_${1}$(printf "%02d" $3)_${2}$(printf "%02d" $4).nc
}
###############################################################################
#                 Clean diags or MMO files after postprocessing               #
#                                                                             #
# $1 : starting date                                                          #
# $2 : expid                                                                  #
# $3 : member                                                                 #
# $4 : starting leadtime                                                      # 
# $5 : end leadtime                                                           #
# $6 : diags/MMO                                                              #
# $7 : list of files extracted from the original tarballs
#                                                                             #
#         Created in May 2012           Author : vguemas@ic3.cat              #
#         Modified: June 2013           isabel.andreu-burillo@ic3.cat         #
###############################################################################

function clean_diagsMMO {
typeset var yyyy0=`echo $1|cut -c1-4`
typeset var mm0=`echo $1|cut -c5-6`
typeset var year1=$(($yyyy0+(10#$mm0+$4-2)/12))
typeset var year2=$(($yyyy0+(10#$mm0+$5-2)/12))
typeset var mon1=$(((10#$mm0+$4-2)%12+1))
typeset var mon2=$(((10#$mm0+$5-2)%12+1))

typeset var listroots
  case $6 in
    'diags' ) listroots="t3d" ;;
    'MMO' ) 
      if [ -z "${7}" ] ; then
        listroots="grid_T grid_U grid_V grid_W icemod"
      else
        listroots=${7}
      fi
     ;;
  esac
typeset var root
typeset var lstfiles
for root in ${listroots[@]} ; do
  rm -f ${root}_$2_$1_fc$3_${year1}$(printf "%02d" $mon1)_${year2}$(printf "%02d" $mon2).nc
done
}
###############################################################################
#                     Vertically averaged salt content                        #
#                                                                             #
# $1 : input grid_T file name                                                 #
# $2 : upper depth of the layer (in meters)                                   #
# $3 : lower depth of the layer (in meters)                                   #
# $4 : output file name  (=> 2D)                                              #
#                                                                             #
#         Created in February 2012      Author : vguemas@ic3.cat              #
###############################################################################

function vertmeansal {
cdo_version=`cdo -V &> ff ; grep Climate ff | cut -d \  -f 5`
rm ff
typeset var ntime=`cdo ntime $1`
typeset var list=""
typeset var jt
for jt in $(seq 1 $ntime); do
  ncks -d time,$((jt-1)) $1 intvertmeansal.nc
  cdfvertmean intvertmeansal.nc vosaline T $2 $3
  ncrename -O -v sovertmean,vertmeansal -d time_counter,time -v time_counter,time vertmean.nc
  mv vertmean.nc outputvertmeansal_$jt.nc
  list=$list" "outputvertmeansal_$jt.nc
  rm -f intvertmeansal.nc
  #test on cdo version: if >1.5.6, remove valid_min/max attributes to avoid values out of that range to be replaced by NaN
  if [[ "$cdo_version" = "`echo -e "$cdo_version\n1.5.6" | sort -V | head -n1`" ]] ; then
   ncatted -O -a valid_max,vertmeansal,d,, outputvertmeansal_$jt.nc outputvertmeansal_$jt.nc 
   ncatted -O -a valid_min,vertmeansal,d,, outputvertmeansal_$jt.nc outputvertmeansal_$jt.nc
  fi
done
cdo cat $list $4
ncks -A -v time $1 $4
rm -f $list 
setminmax $4 vertmeansal
}
###############################################################################
#                    Compute mixed layer heat and salt content                #
#                                                                             #
# $1 : input grid_T file name                                                 #
# $2 : output file name  (=> 2D x-y )                                         #
#                                                                             #
#         Created in February 2012      Author : vguemas@ic3.cat              #
################################################################################

function heat_sal_mxl {
typeset var ntime=`cdo ntime $1`
typeset var list=""
typeset var jt 
typeset var lstvars=`cdo showvar $1`
for jt in $(seq 1 $ntime); do
 ncks -d time,$((jt-1)) $1 intheat_sal_mxl.nc
 if [[ ${lstvars/somxl010} == ${lstvars} ]] ; then
  cdfmxl intheat_sal_mxl.nc mxl.nc
  ncrename -d time_counter,time mxl.nc
  ncks -A mxl.nc intheat_sal_mxl.nc
  rm -f mxl.nc
 fi
 cdfmxlheatc intheat_sal_mxl.nc 
 if [[ $lstvars != ${lstvars/vosaline} ]] ; then  
   cdfmxlsaltc intheat_sal_mxl.nc
   ncks -A mxlsaltc.nc mxlheatc.nc  
   rm -f mxlsaltc.nc
 fi
 mv mxlheatc.nc outputintheat_sal_mxl_$jt.nc
 timevar=`ncdump -h outputintheat_sal_mxl_$jt.nc | grep UNLIMITED | awk '{print $1}'`
 if [[ $timevar == 'time_counter' ]] ; then ncrename -v time_counter,time -d time_counter,time outputintheat_sal_mxl_$jt.nc ; fi
 list=$list" "outputintheat_sal_mxl_$jt.nc
 rm -f intheat_sal_mxl.nc
done
cdo cat $list $2
ncks -A -v time $1 $2
rm -f $list
setminmax $2 somxlheatc
if [[ $lstvars != ${lstvars/vosaline} ]] ; then setminmax $2 somxlsaltc ; fi
}
###############################################################################
#          Pointwise Ocean Heat Content in a specified ocean thickness        #
#          (J/m-2)
#                                                                             #
# $1 : input grid_T file name                                                 #
# $2 : upper depth of the layer (in meters)                                   #
# $3 : lower depth of the layer (in meters)                                   #
# $4 : output file name  (=> 2D x-y )                                         #
#                                                                             #
# Created in June 2012      Author : isabel.andreu-burillo@ic3.cat            #
# May 2014 - Virginie Guemas - Way around the bc that does not work on moore  #
###############################################################################

function ohc_specified_layer {
typeset var ntime=`cdo ntime $1`
typeset var list=""
typeset var jt
ncap2 -v -O -s "heatc_sl=tmask*e3t" mesh_zgr.nc e3t_file.nc
ncrename -d t,time -d z,deptht e3t_file.nc
for jt in $(seq 1 $ntime); do
  cdo seltimestep,$jt $1 intohc_slayer.nc
  ncks -O -v votemper intohc_slayer.nc intmeantem.nc
  ncrename -v votemper,heatc_sl intmeantem.nc #to be commented
  cdo mul intmeantem.nc e3t_file.nc heatc_sl_out.nc
#?  ncks -A -m -v nav_lon,nav_lat $1 heatc_sl_out.nc
  # extract the data between the two given depths --> heatc_sl_top.nc
  ncks -d deptht,$2,$3 heatc_sl_out.nc heatc_sl_top.nc
  #perform the integration of ohc down to that level (main contribution)
  ncap2 -O -s 'heatc_sl=heatc_sl.total($deptht)' heatc_sl_top.nc heatc_sl_top.nc
  # now extract a few levels below, to compute the residual ohc 
  # lower_bnd=`echo "$3 + 200.0" | bc` -> does not work on new moore
  # strip out the .* from $3:
  stripped=`echo ${3/.*}`
  # addition with float returned:
  lower_bnd=`echo $(printf "%f" $(( $stripped + 200)))`  
  ncks -d deptht,$3,$lower_bnd heatc_sl_out.nc heatc_sl_bottom.nc
  # obtain the weight for the extra level containing the 300 m
  # deptht in the gridT files is positive
  # weight = (300.0 - depth_top)/(depth_bottom - depth_top)
  # and add the thickness down to 300 m in the next layer
  ncpdq -a '-deptht' heatc_sl_top.nc heatc_sl_top_invert.nc
  ncks -d deptht,0,0,1 heatc_sl_top_invert.nc level_above.nc
  ncks -d deptht,0,0,1 heatc_sl_bottom.nc level_below.nc
  ## Here, add the residual contribution, before adding it to the main contribution 
  ncrename -v deptht,layerthcknss level_below.nc
  ncrename -v deptht,layerthcknss level_above.nc
  ncbo -A --op_typ=sub -v layerthcknss level_below.nc level_above.nc depth_diff_lay.nc
  ncrename -v layerthcknss,heatc_sl depth_diff_lay.nc
  ncap2 -s "heatc_sl=($3 - layerthcknss)" level_above.nc depth_diff_sublay.nc
  ncbo --op_typ=/ -v heatc_sl depth_diff_sublay.nc depth_diff_lay.nc factor.nc
  ncrename -v heatc_sl,factor factor.nc #to be commented
  ncks -A -v factor factor.nc level_below.nc
  rm -f depth_diff_sublay.nc depth_diff_lay.nc
  ncap2 -O -s "heatc_sl=(factor * heatc_sl)" level_below.nc level_below.nc
  ncwa -O -a deptht level_below.nc level_below.nc
  ncbo --op_typ=+ -v heatc_sl heatc_sl_top.nc level_below.nc total_heatc_sl.nc
  ncap2 -s "heatc_sl=1020.0*4000*heatc_sl" total_heatc_sl.nc heatc_sl_$jt.nc
  list=$list" "heatc_sl_$jt.nc
  rm -f depth_diff_lay.nc depth_diff_sublay.nc
  rm -f heatc_sl_out.nc heatc_sl_top.nc heatc_sl_top_invert.nc heatc_sl_bottom.nc 
  rm -f level_above.nc level_below.nc
  rm -f intohc_slayer.nc intmeantem.nc vertmean.nc total_heatc_sl.nc
  rm -f factor.nc
done
cdo cat $list $4
ncks -A -v time $1 $4
rm -f $list 
rm -f e3t_file.nc
setminmax $4 heatc_sl
}
###############################################################################
#                      Compute the MOC for oceanic basins                     #
#                                                                             #
# $1 : input grid_V file name                                                 #
# $2 : output file name (=> 2D, depth-y)                                      #
#                                                                             #
#         Created in March 2012         Author : vguemas@ic3.cat              #
###############################################################################

function moc {
typeset var ntime=`cdo ntime $1`
typeset var list=""
typeset var jt
for jt in $(seq 1 $ntime); do
 cdo seltimestep,$jt $1 intmoc.nc
 cdfmoc intmoc.nc
 ncwa -O -a x moc.nc outmoc_$jt.nc
 ncks -O -x -v nav_lon,nav_lat outmoc_$jt.nc outmoc_$jt.nc
 timevar=`ncdump -h outmoc_$jt.nc | grep UNLIMITED | awk '{print $1}'`
 if [[ $timevar == 'time_counter' ]] ; then ncrename -v time_counter,time -d time_counter,time outmoc_$jt.nc ; fi
 list=$list" "outmoc_$jt.nc
 rm -f intmoc.nc moc.nc
done
cdo cat $list $2
lstdims=`ncdump -h $2 | awk /dimensions:/,/variables:/ | grep -v dimensions: | grep -v variables: | awk '{print $1}'`
if [[ ${lstdims/gsize} != ${lstdims} ]] ; then
  ncrename -d gsize,y $2
fi
ncks -A -v nav_lon,nav_lat $1 $2
ncks -A -v time $1 $2
rm -f $list
}
###############################################################################
#                                                                             #
#   Compute the intensity of convection in the four main convection sites     #
#                                                                             #
# $1 : input oce file name containing somxl010                                #
# $2 : input grid                                                             #
# $3 : output file name (=> index)                                            #
#                                                                             #
#         Created in October 2013         Author : vguemas@ic3.cat            #
###############################################################################

function convection {
case $2 in
  'Ec2.3_O1L42'|'Ec3.0_O1L46'|'N3.2_O1L42'|'N3.3_O1L46'|'nemovar_O1L42') 
    A1=225;A2=245;A3=215;A4=255;
    B1=245;B2=290;B3=215;B4=245;
    C1=260;C2=310;C3=245;C4=291;
    D1=225;D2=280;D3=1;D4=50;;

  'Ec3.0_O25L46'|'Ec3.0_O25L75'|'glorys2v1_O25L75')
    stop"Option convection not available yet for this configuration"
  ;;
esac

cdo fldmax -selindexbox,${A1},${A2},${A3},${A4} $1 Labrador.nc
ncrename -v somxl010,Labrador Labrador.nc
ncks -v Labrador Labrador.nc convection.nc
rm -f Labrador.nc 

cdo fldmax -selindexbox,${B1},${B2},${B3},${B4} $1 Irminger.nc
ncrename -v somxl010,Irminger Irminger.nc
ncks -A -v Irminger Irminger.nc convection.nc
rm -f Irminger.nc 

cdo fldmax -selindexbox,${C1},${C2},${C3},${C4} $1 GIN.nc
ncrename -v somxl010,GIN GIN.nc
ncks -A -v GIN GIN.nc convection.nc
rm -f GIN.nc 

cdo fldmax -selindexbox,${D1},${D2},${D3},${D4} $1 Wedell.nc
ncrename -v somxl010,Wedell Wedell.nc
ncks -A -v Wedell Wedell.nc convection.nc
rm -f Wedell.nc 

mv convection.nc $3
}
###############################################################################
#                                                                             #
#                     Compute the barotropic stream function                  #
#                                                                             #
# $1 : input grid_U file name                                                 #
# $2 : input grid_V file name                                                 #
# $3 : output file name without nc extension  (=> 2D x-y)                     #
#                                                                             #
#         Created in March 2012         Author : vguemas@ic3.cat              #
###############################################################################

function psi {
typeset var ntime=`cdo ntime $1`
typeset var list=""
typeset var jt
for jt in $(seq 1 $ntime); do
  cdo seltimestep,$jt $1 intU.nc
  cdo seltimestep,$jt $2 intV.nc
  cdfpsi intU.nc intV.nc
  mv psi.nc psi_U.nc
  cdfpsi intU.nc intV.nc V 
  mv psi.nc psi_V.nc
  ncea psi_U.nc psi_V.nc psi_${jt}.nc
  timevar=`ncdump -h psi_$jt.nc | grep UNLIMITED | awk '{print $1}'`
  if [[ $timevar == 'time_counter' ]] ; then ncrename -v time_counter,time -d time_counter,time psi_$jt.nc ; fi
  list=$list" "psi_$jt.nc
  rm -f intU.nc intV.nc psi_U.nc psi_V.nc
done
cdo cat $list ${3}
ncks -A -v time $1 ${3}
rm -f $list
}
###############################################################################
#                                                                             #
#        Compute the intensity of the subtropical and subpolar gyres          #
#                                                                             #
# $1 : input psi file name                                                    #
# $2 : input grid                                                             #
# $3 : output file name   ( => index )                                        #
#                                                                             #
#         Created in October 2013         Author : vguemas@ic3.cat            #
###############################################################################

function gyres {
case $2 in
  'Ec2.3_O1L42'|'Ec3.0_O1L46'|'N3.2_O1L42'|'N3.3_O1L46'|'nemovar_O1L42') 
    A1=230;A2=275;A3=215;A4=245;
    B1=70;B2=145;B3=195;B4=235;
    C1=45;C2=175;C3=165;C4=220;
    D1=195;D2=275;D3=175;D4=225; 
    E1=70;E2=205;E3=120;E4=145;
    F1=235;F2=300;F3=120;F4=145;
    G1=320;G2=30;G3=110;G4=180;
    H1=1;H2=361;H3=1;H4=65;;

  'Ec3.0_O25L46'|'Ec3.0_O25L75'|'glorys2v1_O25L75')
    stop"Option gyres not available yet for this configuration"
  ;;
esac

cdo fldmin -selindexbox,${A1},${A2},${A3},${A4} $1  subpolar_NAtl.nc
ncrename -v sobarstf,subpolNAtl subpolar_NAtl.nc
cdo mulc,-1 subpolar_NAtl.nc gyres.nc
rm -f subpolar_NAtl.nc

cdo fldmin -selindexbox,${B1},${B2},${B3},${B4} $1 subpolar_NPac.nc
ncrename -v sobarstf,subpolNPac subpolar_NPac.nc
cdo mulc,-1 subpolar_NPac.nc tmp.nc
ncks -A tmp.nc gyres.nc
rm -f subpolar_NPac.nc tmp.nc

cdo fldmax -selindexbox,${C1},${C2},${C3},${C4} $1 subtrop_NPac.nc
ncrename -v sobarstf,subtropNPac subtrop_NPac.nc
ncks -A subtrop_NPac.nc gyres.nc
rm -f subtrop_NPac.nc

cdo fldmax -selindexbox,${E1},${E2},${E3},${E4} $1 subtrop_SPac.nc
ncrename -v sobarstf,subtropSPac subtrop_SPac.nc
ncks -A subtrop_SPac.nc gyres.nc
rm -f subtrop_SPac.nc

cdo fldmax -selindexbox,${D1},${D2},${D3},${D4} $1 subtrop_NAtl.nc
ncrename -v sobarstf,subtropNAtl subtrop_NAtl.nc
ncks -A subtrop_NAtl.nc gyres.nc
rm -f subtrop_NAtl.nc

cdo fldmax -selindexbox,${F1},${F2},${F3},${F4} $1 subtrop_SAtl.nc
ncrename -v sobarstf,subtropSAtl subtrop_SAtl.nc
ncks -A subtrop_SAtl.nc gyres.nc
rm -f subtrop_SAtl.nc

cdo fldmax -selindexbox,${G1},${G2},${G3},${G4} $1 subtrop_Ind.nc
ncrename -v sobarstf,subtropInd subtrop_Ind.nc
ncks -A subtrop_Ind.nc gyres.nc
rm -f subtrop_Ind.nc

cdo fldmax -selindexbox,${H1},${H2},${H3},${H4} $1 ACC.nc
ncrename -v sobarstf,ACC ACC.nc
ncks -A ACC.nc gyres.nc
rm -f ACC.nc

mv gyres.nc $3

}
###############################################################################
#                                                                             #
#   Compute an Atlantic MOC index by averaging the meridional overturning     #
#                in a latitude band between 1km and 2km                       #
#        or any other index averaging the meridional overturning in           #
#                   a given basin and a given domain                          #
#                                                                             #
# $1 : input moc file name                                                    #
# $2 : latitude min                                                           #
# $3 : latitude max                                                           #
# $4 : output file name      ( => index )                                     #
# $5 : depth min (default : 1km)                                              #
# $6 : depth max (default : 2km)                                              #
# $7 : basin (default :  zomsfatl)                                            # 
#                                                                             #
#         Created in March 2012         Author : vguemas@ic3.cat              #
###############################################################################

function area_moc {
if [ -z "$5" ] ; then
  typeset var depmin=-1000.0
else
  typeset var depmin=-$5
fi
if [ -z "$6" ] ; then
  typeset var depmax=-2000.0
else
  typeset var depmax=-$6
fi
if [ -z "$7" ] ; then
  typeset var basin=zomsfatl
else
  typeset var basin=$7
fi
lstdims=`ncdump -h $1 | awk /dimensions:/,/variables:/ | grep -v dimensions: | grep -v variables: | awk '{print $1}'`
if [[ ${lstdims/x} != ${lstdims} ]] ; then
  ncwa -O -a x $1 tmpmoc.nc
else
  cp $1 tmpmoc.nc
fi
ncrename -O -d y,lat -v nav_lat,lat tmpmoc.nc tmpmoc.nc
ncks -O -v $basin,time,depthw,lat  tmpmoc.nc  tmpmoc.nc 
ncks -d lat,$2,$3 -d depthw,${depmax},${depmin} tmpmoc.nc area_moc.nc
cdo vertmean area_moc.nc area_ave_moc.nc
ncap -O -s "coslat[lat]=cos(lat[lat]*3.141592657/180.0)" area_ave_moc.nc area_ave_moc2.nc
ncwa -w coslat -a lat area_ave_moc2.nc area_ave_moc3.nc
ncks -O -v $basin,time area_ave_moc3.nc $4
rm -f tmpmoc.nc area_moc.nc area_ave_moc2.nc area_ave_moc3.nc
if [[ $4 != area_ave_moc.nc ]] ; then
  rm -f area_ave_moc.nc
fi
}
###############################################################################
#                                                                             #
#   Compute an Atlantic MOC index by finding the maximum  of the annual       #
#        mean meridional overturning in a latitude / depth region             #
#                                                                             #
# $1 : input moc file name                                                    #
# $2 : latitude min                                                           #
# $3 : latitude max                                                           #
# $4 : depth mean                                                             #
# $5 : depth max                                                              #
# $6 : output file name        ( => index )                                   #
#                                                                             #
#         Created in March 2012         Author : vguemas@ic3.cat              #
###############################################################################

function max_moc {
if [ ! -f $6 ] ; then
 ncecat -h $1 tmpmoc1.nc
 lstdims=`ncdump -h tmpmoc1.nc | awk /dimensions:/,/variables:/ | grep -v dimensions: | grep -v variables: | awk '{print $1}'`
 if [[ ${lstdims/x} != ${lstdims} ]] ; then
   ncwa -O -a x tmpmoc1.nc tmpmoc1.nc
 fi
 ncrename -d record,x tmpmoc1.nc
 ncpdq -O -h -a time,x tmpmoc1.nc tmpmoc1.nc
 ncpdq -O -h -a depthw,x tmpmoc1.nc tmpmoc1.nc
 ncpdq -O -h -a y,x tmpmoc1.nc tmpmoc1.nc
 cdo yearmean tmpmoc1.nc tmpmoc.nc 
 typeset var ntime=`cdo ntime tmpmoc.nc`
 typeset var list=""
 for jt in $(seq 1 $ntime) ; do
   cdo seltimestep,$jt tmpmoc.nc tmpmoc2.nc
   cdfmaxmoc tmpmoc2.nc atl $2 $3 $4 $5
   mv maxmoc.nc maxmoc_$jt.nc
   timevar=`ncdump -h maxmoc_$jt.nc | grep UNLIMITED | awk '{print $1}'`
   if [[ $timevar == 'time_counter' ]] ; then ncrename -v time_counter,time -d time_counter,time maxmoc_$jt.nc ; fi
   list=${list}" "maxmoc_$jt.nc
   rm -f tmpmoc2.nc
 done
 cdo cat $list $6
 ncks -A -v time tmpmoc.nc $6
 rm -f $list tmpmoc.nc tmpmoc1.nc
fi
}
###############################################################################
#                                                                             #
#        Compute the sia extent, area and volume in both hemispheres          #
#                                                                             #
# $1 : input ice file name                                                    #
# $2 : output file name     ( => index )                                      #
#                                                                             #
#         Created in April 2012         Author : vguemas@ic3.cat              #
###############################################################################

function siasiesiv {
cp ${CON_FILES}/ice_template.nc toto_N.nc
cp ${CON_FILES}/ice_template.nc toto_S.nc

case ${NEMOVERSION} in
  'Ec3.0_O1L46'|'Ec3.0_O25L46'|'Ec3.0_O25L75') for var in `cdo showvar $1 | head -1`
do
[[ $var = "ice_pres" || $var = "iiceconc" ]] && ncrename -v $var,ileadfra $1 
done;;
#'Ec3.0_O1L46'|'Ec3.0_O25L46') ncrename -v ice_pres,ileadfra $1 ;;
#'Ec3.0_O1L46'|'Ec3.0_O25L46') ncrename -v iiceconc,ileadfra $1 ;;
esac
  
typeset var ntime=`cdo ntime $1`
typeset var list1=""
typeset var list2=""
typeset var jt
for jt in $(seq 1 $ntime) ; do
  cdo seltimestep,$jt $1 tmpice.nc
  cdficediags tmpice.nc>ice.txt
  for d in N S;do
    ncdump toto_${d}.nc  > ice_template.cdl
    sia=`grep ${d}Area ice.txt |awk '{print $4}'`
    sie=`grep ${d}Exnsidc ice.txt|awk '{print $4}'`
    siv=`grep ${d}Volume ice.txt|awk '{print $4}'`
    sed -e "s/sia =.*/sia = $sia ;/" ice_template.cdl > ice_template2.cdl
    sed -e "s/sie =.*/sie = $sie ;/" ice_template2.cdl > ice_template3.cdl
    sed -e "s/siv =.*/siv = $siv ;/" ice_template3.cdl > ice_template.cdl
    ncgen -o ice_${d}_${jt}.nc ice_template.cdl
    rm -f ice_template.cdl ice_template2.cdl ice_template3.cdl
  done
  list1=$list1" "ice_N_${jt}.nc
  list2=$list2" "ice_S_${jt}.nc
  rm -f ice.txt tmpice.nc icediags.nc
done
cdo cat $list1 ice_N_${2}
cdo cat $list2 ice_S_${2}
ncks -A -v time $1 ice_N_${2}
ncks -A -v time $1 ice_S_${2}
rm -f $list1 $list2 toto_N.nc toto_S.nc
setminmax ice_N_${2} sia sie siv
setminmax ice_S_${2} sia sie siv
}
###############################################################################
#                                                                             #
#                    Compute the total ocean heat extent                      #
#                                                                             #
# $1 : input temperature file name                                            #
# $2 : output file name  ( => 2D x-y )                                        #
# $3 : basin (NAtl, NPac, TAtl, TPac, TInd, Anta, Arct, Glob) Default : Glob  #
#      $4 = 0 if $3 = Glob                                                    #
# $4 : mixed layer (1=only, 0=included, -1=without) Default : 0               # 
# $5 : upper level of the layer (optional) Default : top                      #
# $6 : lower level of the layer (optional) Default : bottom                   #
#                                                                             #
#         Created in May 2012           Author : vguemas@ic3.cat              #
###############################################################################

function ohc {
cp ${CON_FILES}/depth.${NEMOVERSION}.txt depth.txt
#
# Input arguments
#
if [ -z "$3" ] ; then
  typeset var basin='Glob'
else
  typeset var basin=$3
fi
if [ -z "$4" ] ; then
  typeset var mxl=0
else
  typeset var mxl=$4
fi
if [ -z "$5" ] ; then 
  typeset var up=1 
else
  typeset var up=$5
fi
if [ -z "$6" ] ; then
  typeset var down=`cat depth.txt | wc -l`  
else
  typeset var down=$6
fi

if [[ ${up} -eq 1 ]] ; then
  typeset var depmin=0
else
  typeset var depmin=`cat depth.txt |head -${up} |tail -1 | awk '{print $2}' | awk '{printf "%.0f",$1}'`
fi
typeset var depmax=`cat depth.txt |head -${down} |tail -1 | awk '{print $2}' | awk '{printf "%.0f",$1}'`

cp ${CON_FILES}/heatc_template.nc template_heatc.nc
ncdump template_heatc.nc > template_heatc.cdl
#exec=/home/huidu/download/CDFTOOLS_2.1/cdfheatc-cfu
exec=/home/vguemas/CDFTOOLS_2.1/cdfheatc-cfu
#
#  Define some parameters
#
typeset var para
typeset var output
typeset var nlev=`cat depth.txt | wc -l`
if [[ ! -z "$depmin" && ! -z "$depmax" ]] ; then
  if [[ $depmin != 0 || ${down} != ${nlev} && ${down} != 0 ]] ; then
    output=${depmin}-${depmax}'_'
  fi
fi

case $basin in
 'NAtl') para="atl $mxl 0 0 10 65"; output='NAtl_10N65N_'${output} ;;
 'TAtl') para="atl $mxl 0 0 -30 30" ; output='TAtl_30S30N_'${output} ;;
 'NPac') para="pac $mxl 0 0 10 70" ; output='NPac_10N70N_'${output} ;;
 'TPac') para="pac $mxl 0 0 -30 30" ; output='TPac_30S30N_'${output} ;;
 'Arct') para="atl $mxl 0 0 65 90" ; output='Arc_65N90N_'${output} ;;
 'Anta') para="all $mxl 0 0 -90 -60" ; output='Ant_90S60S_'${output} ;;
 'TInd') para="ind $mxl 0 0 -30 30" ; output='TInd_30S30N_'${output} ;;
 'Glob') para="all $mxl 0 0 0 0" ;;
esac

case $mxl in
 1) output='mxl_'${output} ;;
 -1) output='nonmxl_'${output} ;; 
esac
#
# Compute ohc
#
typeset var lstvars=`cdo showvar $1`
typeset var ntime=`cdo ntime $1`
typeset var list=""
typeset var jt
for jt in $(seq 1 $ntime) ; do
 cdo seltimestep,$jt $1 tmpohc.nc
 lstdims=`ncdump -h tmpohc.nc | awk /dimensions:/,/variables:/ | grep -v dimensions: | grep -v variables: | awk '{print $1}'`
 if [[ ${lstdims/x_2} != ${lstdims} ]] ; then
   if [[ ${lstdims/x} != ${lstdims} ]] ; then
     ncwa -O -a x tmpohc.nc tmpohc.nc
   fi
   ncrename -d x_2,x tmpohc.nc
 fi
 if [[ ${lstvars/somxl010} != ${lstvars} ]] ; then
  ncks -v somxl010 tmpohc.nc mxl.nc
 else
  cdfmxl tmpohc.nc mxl.nc
 fi
 $exec tmpohc.nc $para $up $down > tmp.log
 cat tmp.log
 thc=`cat tmp.log | grep "Total Heat content :" | awk '{print $5}'`;
 uhc=`cat tmp.log | grep "Total Heat content/volume" | awk '{print $5}'`;
 sed -e "s/thc =.*/thc = $thc ;/" template_heatc.cdl > template_heatc2.cdl
 sed -e "s/uhc =.*/uhc = $uhc ;/" template_heatc2.cdl > template_heatc.cdl
 ncgen -o heatc_${jt}.nc template_heatc.cdl
 rm -f template_heatc2.cdl tmpohc.nc mxl.nc tmp.log
 list=$list" "heatc_${jt}.nc
done
cdo cat $list ${output}$2
ncks -h -A -v time $1 ${output}$2
rm -f $list template_heatc.nc template_heatc.cdl depth.txt
setminmax ${output}$2 thc uhc
}
###############################################################################
#                                                                             #
#                    Cut a meridional or zonal section                        #
#                                                                             #
#                                                                             #
# $1 : input file                                                             #
# $2 : input var                                                              #
# $3 : Z/M   (zonal / meridional section)                                     #
# $4 : lat/lon                                                                #
# $5 : output file  ( => 2D )                                                 #
#                                                                             #
#         Created in September 2012     Author : vguemas@ic3.cat              #
#                                                                             #
###############################################################################

function cutsection {
  typeset var ntime=`cdo ntime $1`
  typeset var nx=`ncdump -h $1|grep 'x = '|head -1|cut -f3 -d" "`
  typeset var ny=`ncdump -h $1|grep 'y = '|head -1|cut -f3 -d" "`
  typeset var nz=`ncdump -h $1|grep 'depth'|head -1|cut -f3 -d" "`
cat>section.R<<EOF1
  library(ncdf) 
  # Get latitudes, longitudes, depth, mask
  fncm=open.ncdf('mesh_hgr.nc')
  lon=get.var.ncdf(fncm,'nav_lon')
  lat=get.var.ncdf(fncm,'nav_lat')
  close.ncdf(fncm)
  fncm=open.ncdf('mesh_zgr.nc')
  depth=get.var.ncdf(fncm,'nav_lev')
  close.ncdf(fncm)
  fncm=open.ncdf('mask.nc')
  mask=get.var.ncdf(fncm,switch('$2','vozocrtx'='umask','vomecrty'='vmask','tmask'))
  close.ncdf(fncm)
  # Latitude / longitude of the zonal / meridional section
  exactpos=$4
  if ('$3'=='M') {
    while (exactpos<min(lon)) {exactpos=exactpos+360}
    while (exactpos>max(lon)) {exactpos=exactpos-360}
  }
  # Collect the indexes defining the section
  listi=array(dim=switch('$3','Z'=$nx-2,'M'=$ny-1))
  listj=array(dim=switch('$3','Z'=$nx-2,'M'=$ny-1))
  for (jpt in 1:length(listi)) {
    vect=switch('$3','Z'=lat[jpt,],'M'=lon[,jpt+1])
    if (min(abs(vect-exactpos))<(2*360./$nx)) {
      pos=sort(abs(vect-exactpos),index.return=T)\$ix[1]
      listi[jpt]=switch('$3','Z'=jpt+1,'M'=pos)
      listj[jpt]=switch('$3','Z'=pos,'M'=jpt)
    }
  }
  listi=listi[is.na(listi)==F]
  listj=listj[is.na(listj)==F]
  print(listi)
  print(listj)
  # Select variable at those indexes
  fnc1=open.ncdf('$1')
  varout=array(dim=c(length(listi),$nz,$ntime))
  for (jt in 1:$ntime) {
    varin=get.var.ncdf(fnc1,'$2',start=c(1,1,1,jt),count=c($nx,$ny,$nz,1))
    varin[which(mask<0.5)]=1e20
    for (jpt in 1:length(listi)) {
      varout[jpt,,jt]=varin[listi[jpt],listj[jpt],]
    }  
  }
  close.ncdf(fnc1)
  # Write the output
  wtime=dim.def.ncdf("time","",seq(1,$ntime),unlim=TRUE) 
  dimout=array(dim=length(listi))
  for (jpt in 1:length(listi)) {
    dimout[jpt]=switch('$3','Z'=lon[listi[jpt],listj[jpt]],'M'=lat[listi[jpt],listj[jpt]])
  }
  wsec=switch('$3','Z'=dim.def.ncdf("lon","",dimout),'M'=dim.def.ncdf("lat","",dimout))
  wdep=dim.def.ncdf("deptht","",depth)
  wvar=var.def.ncdf("$2","",list(wsec,wdep,wtime),1e20)
  fnc2=create.ncdf('$5',wvar)
  put.var.ncdf(fnc2,wvar,varout)
  close.ncdf(fnc2)
EOF1
R CMD BATCH section.R
ncks -h -A -v time $1 $5
}
###############################################################################
#                                                                             #
#   3-dimensional conservative interpolation to the regular atmospheric grid  #
#                                                                             #
# $1 : input file                                                             #
# $2 : input var                                                              #
# $3 : output file         ( => 3D )                                          #
#                                                                             #
#         Created in November 2012      Author : vguemas@ic3.cat              #
#                                                                             #
###############################################################################

function interp3d {
  typeset var nz=`ncdump -h $1|grep 'deptht'|head -1|cut -f3 -d" "`
  ln -sf  /cfu/pub/scripts/interpolation/scrip_use scrip_use
  for lev in $(seq 1 $nz) ; do
    ncks -d deptht,$((lev-1)) -v $2 $1 tmp_${lev}.nc
    ncwa -O -h -a deptht tmp_${lev}.nc tmp_${lev}.nc
    ln -sf /cfu/autosubmit/con_files/weigths/${NEMOVERSION}/rmp_${NEMOVERSION}_to_*_lev${lev}.nc rmp_${NEMOVERSION}_to_regular_lev${lev}.nc
    cat > scrip_use_in <<EOF
&remap_inputs
    remap_wgt   = 'rmp_${NEMOVERSION}_to_regular_lev${lev}.nc'
    infile      = 'tmp_${lev}.nc'
    invertlat   = FALSE
    var         = '$2'
    fromregular = FALSE
    outfile     = 'tmp_$(printf "%02d" $lev).nc'
/
EOF
    ./scrip_use
    ncecat -O -h tmp_$(printf "%02d" $lev).nc tmp_$(printf "%02d" $lev).nc
    ncrename -h -O -d record,deptht tmp_$(printf "%02d" $lev).nc tmp_$(printf "%02d" $lev).nc 
  done
  ncrcat -n $nz,2,1 -v $2 tmp_01.nc tmpout.nc
  ncpdq -O -h -a time,deptht tmpout.nc tmp2out.nc 
  ncks -A -v deptht $1 tmp2out.nc
  cdo setgrid,t106grid tmp2out.nc tmp3out.nc
  orca_res=`echo $NEMOVERSION | cut -c 7-9`
  if ((${orca_res} \= "O25")); then
   cdo invertlatdata tmp3out.nc $3
  else
    cp tmp3out.nc $3
  fi
  rm -f tmp_* tmpout.nc tmp2out.nc tmp3out.nc rmp_${NEMOVERSION}_to_regular_lev* scrip_use scrip_use_in
}

###############################################################################
#                                                                             #
#                    Set minimum and maximum attributes                       #
#                                                                             #
# $1 : input file                                                             #
# $2 : input vars provided as a list                                          #
#                                                                             #
#         Created in May 2012           Author : vguemas@ic3.cat              #
###############################################################################

function setminmax {
  cdo fldmax -timmax $1 intmax.nc
  cdo fldmin -timmin $1 intmin.nc
  typeset var listvars=""
  typeset args=("$@")
  typeset var ind
  for ((ind=1;ind<$#;ind++)) ; do
    listvars=$listvars" "${args[$ind]}
  done
  typeset name
  for name in ${listvars[@]} ; do
    typeset var var_max=`ncdump -v $name intmax.nc | tail -2 | head -1 | awk '{print $1}'`
    ncatted -h -a valid_max,$name,m,f,$var_max $1
    typeset var var_min=`ncdump -v $name intmin.nc | tail -2 | head -1 | awk '{print $1}'`
    ncatted -h -a valid_min,$name,m,f,$var_min $1
  done
  rm -f intmax.nc intmin.nc
}

###############################################################################
#                                                                             #
#              Concatenate with a weight of the results obtained              #
#                             on partial years                                #
#                                                                             #
# $1 : input file 1                                                           #
# $2 : input file 2                                                           #
# $3 : ending month file 1                                                    # 
# $4 : output file                                                            #
#                                                                             #
#         Created in May 2012           Author : vguemas@ic3.cat              #
###############################################################################

function concat {
  typeset var wgt=$((10#${3}%12))
  typeset var ntimes=`cdo ntime $1`
  ncks -O -d time,$(($ntimes-1)) $1 part1.nc
  ncks -O -d time,0 $2 part2.nc
  cdo mulc,$wgt part1.nc part1_bis.nc
  cdo mulc,$((12-$wgt)) part2.nc part2_bis.nc
  cdo add part1_bis.nc part2_bis.nc tmp2_bis.nc
  cdo divc,12 tmp2_bis.nc tmp2.nc
  lstdims=`ncdump -h tmp2.nc | awk /dimensions:/,/variables:/ | grep -v dimensions: | grep -v variables: | awk '{print $1}'`
  if [[ ${lstdims/lev} != ${lstdims} ]] ; then
    ncrename -d lev,ensemble tmp2.nc 
  else
    ncecat -O -h tmp2.nc tmp2.nc
    ncrename -h -O -d record,ensemble tmp2.nc tmp2.nc
    ncpdq -O -h -a time,ensemble tmp2.nc tmp2.nc
  fi
  rm -f part1.nc part2.nc part1_bis.nc part2_bis.nc tmp2_bis.nc
  ncks -O -d time,0,$((ntimes-2+10#${3}/12)) $1 tmp1.nc
  ncks -O -d time,1, $2 tmp3.nc
  ncrcat -O tmp1.nc tmp2.nc tmp3.nc $4
  rm -f tmp1.nc tmp2.nc tmp3.nc
}

###############################################################################
#                                                                             #
#                Gather the members in a single netcdf file                   #
#                                                                             #
# $1 : prefix netcdf file name for all the members                            #
# $2 : suffix netcdf file name for all the members                            #
# $3 : output file name                                                       #
#                                                                             #
#         Created in May 2012           Author : vguemas@ic3.cat              #
###############################################################################        

function gather_memb {
  typeset var lstfiles1=`ls ${1}*${2}`
  typeset var file
  typeset var ind=1
  typeset var lstfiles2=""
  for file in ${lstfiles1[@]} ; do
    typeset var lstvars=`cdo showvar $file`
    typeset var lstexclude="nav_lon nav_lat"
    for varex in $lstexclude ; do
      if [[ ${lstvars/${varex}/} != ${lstvars} ]] ; then
        ncks -O -x -v ${varex} $file $file
      fi
    done
    if [[ ! -z "`ncdump -h ${file} | grep ensemble | head -1 | awk '{print $3}' | grep [0-9]`" ]] ; then
      typeset var nrec=`ncdump -h ${file} | grep ensemble | head -1 | awk '{print $3}' | grep [0-9]`
      typeset var jrec
      for jrec in $(seq 1 $nrec) ; do
        ncks -d ensemble,$((jrec-1)) $file tmp_record_${ind}.nc
        ncwa -O -a ensemble tmp_record_${ind}.nc tmp_record_${ind}.nc
        lstfiles2=${lstfiles2}" "tmp_record_${ind}.nc
        ind=$((ind+1))
      done
    else 
      lstfiles2=${lstfiles2}" "$file
    fi
  done
  ncecat -h ${lstfiles2} ${3}
  rm -f tmp_record_*.nc
  ncrename -h -O -d record,ensemble ${3} ${3}
  ncpdq -O -h -a time,ensemble ${3} ${3}
}

###############################################################################
#          Choose vertical level in ocean, or vertically average between 
#           2 or more  ocean levels                                           #
#                                                                             #
# $1 : input file name                                                        #
# $2:  input variable name (votemper or vosaline)                             #
# $3 : upper depth of the layer (in level number)                             #         
# $4 : lower depth of the layer (in level number. If same as the              # 
#      upper layer, then the script just selects this layer. If its           #
#      different, a vertical integral (weighted) between upper and lower      #
#      level is computed                                                      #
# $5 : output file name  (=> 2D field  )                                      #
#                                                                             #
# Created in  February 2012      Author : vguemas@ic3.cat                     #
# Modified (more generic, 
# i.e. for any input var) in December 2014                                    #
#                                Author : eleftheria.exarchou@ic3.cat         #
###############################################################################

function vertmeanvar {
 typeset var ntime=`cdo ntime $1`
 typeset var list=""
 typeset var jt
for jt in $(seq 1 $ntime); do
  ncks -d time,$((jt-1)) $1  a1
 # The oras4 data do not have gdepth data in their mask, but only gdept_0, so: 
  if    grep  -q  nemovar_s4  $1   ; then 
  l1=$(( $3 % 6 )) 
  l2=$(( $4 % 6 )) 
  ll1=$(( $3 / 6 + 1 )) 
  ll2=$(( $4 / 6 + 1 )) 
  lev1=`echo $(cdo  output  -selvar,gdept_0   mesh_zgr.nc   |  sed -n ${ll1}p  |  awk '{ print $'$l1' }')`
  lev2=`echo $(cdo  output  -selvar,gdept_0   mesh_zgr.nc   |  sed -n ${ll2}p  |  awk '{ print $'$l2' }')`
  else
  l1=$(($3+1))
  l2=$(($4+1))
  lev1=`echo $(cdo info  -seltimestep,1    -selvar,gdept   mesh_zgr.nc   |  sed -n ${l1}p  |  awk '{ print $10 }')`
  lev2=`echo $(cdo info  -seltimestep,1    -selvar,gdept   mesh_zgr.nc   |  sed -n ${l2}p  |  awk '{ print $10 }')`
  fi 
  cdfvertmean a1  $2 T $lev1 $lev2
  rm -f a1 
  ncrename -O -v sovertmean,vertmean -d time_counter,time -v time_counter,time vertmean.nc
  mv vertmean.nc outputvertmean_$jt.nc
  list=$list" "outputvertmean_$jt.nc
  rm -f intvertmean.nc a? 
#  #test on cdo version: if >1.5.6, remove valid_min/max attributes to avoid values out of that range to be replaced by NaN
   ncatted -O -a valid_max,vertmean,d,,                             outputvertmean_$jt.nc  outputvertmean_$jt.nc 
   ncatted -O -a valid_min,vertmean,d,,                             outputvertmean_$jt.nc  outputvertmean_$jt.nc
   ncatted -O -a standard_name,time,a,c,time                        outputvertmean_$jt.nc  outputvertmean_$jt.nc
   ncatted -O -a units,time,o,c,'seconds since 1993-05-01 00:00:00' outputvertmean_$jt.nc  outputvertmean_$jt.nc
   ncatted -O -a long_name,time,a,c,'Time axis'                     outputvertmean_$jt.nc  outputvertmean_$jt.nc
done
cdo cat $list $5
ncks -A -v time $1 $5
rm -f $list   
setminmax $5 vertmean
ncrename -v vertmean,$2 $5 
#typeset var level=`echo $(cdo info  -selvar,$2  -setctomiss,0  $5   |  sed -n 2p | awk '{ print $7 }')`
#typeset var lev=`echo $(cdo info  -seltimestep,1  -selvar,$2  -setctomiss,0  $1   | grep $level  | awk '{ print $1 }')`
lev=$3
echo $lev
cp $5  tmp_${lev}.nc
## Here we interpolate horizontally onto a regular grid
    ln -sf /cfu/autosubmit/con_files/weigths/${NEMOVERSION}/rmp_${NEMOVERSION}_to_*_lev${lev}.nc rmp_${NEMOVERSION}_to_regular_lev${lev}.nc
    ln -sf  /cfu/pub/scripts/interpolation/scrip_use scrip_use
    cat > scrip_use_in <<EOF
&remap_inputs
    remap_wgt   = 'rmp_${NEMOVERSION}_to_regular_lev${lev}.nc'
    infile      = 'tmp_${lev}.nc'
    invertlat   = FALSE
    var         = '$2'
    fromregular = FALSE
    outfile     = 'tmp_${lev}_int.nc'
/
EOF
    ./scrip_use
   cdo setgrid,t106grid tmp_${lev}_int.nc tmpout1.nc
   cdo invertlatdata tmpout1.nc  $5
   ncatted -O -a long_name,${2},o,c,'vertical average of temp/sal between '${lev1}'m and '${lev2}'m'  $5
   rm -f tmp_* tmpout*  rmp_${NEMOVERSION}_to_regular_lev* scrip_use scrip_use_in
}