#!/bin/bash
#
# computes and plots some simple ensemble statistics for the HIPPO_r0xxi1850p1f12xx ensemble
# 
# to be run as a cronjob on climcal3.
#
# ralf.hand@giub.unibe.ch 2020-03-25

set -e 

dir=/mnt/climstor/ERC_PALAEO

setlist="HIPPO_set_1850_1"

cdo selvar,tsurfgm HIPPO_r001i1850p1f1200/means/HIPPO_r001i1850p1f1200_1850_mon.special.nc gridref.nc
cdo griddes gridref.nc >> gridref_timeseries.dat
gridref_timeseries=$(readlink -f gridref_timeseries.dat)

oldpwd=$(pwd)

for set in ${setlist}; do
   
   ensstatdir=${dir}/${set}/enststat
   
   echo ${ensemble}

   if [ ! -e ${ensstatdir} ]; then mkdir -p ${ensstatdir}; fi
   
   tmp=$(mktemp -d /scratch3/ralf/spielwiese/tmp.ensstat.XXX)
   cd ${tmp}
   
   for year in $(seq 1850 2009); do
   
      echo "computing enstat for year ${year}"
  
      set +e

      filelist=$(ls ${dir}/${set}/HIPPO_r*i*p*f*/means/HIPPO_r*i*p*f*_${year}_mon.special.nc 2>/dev/null)

      set -e

      echo ${filelist}

      ensstat_filename=${set}_ENSSTAT_TS.special.nc

      nfiles=$(echo ${filelist} | wc -w)
      echo "number of files for year ${year} ${nfiles}"
   
      for file in ${filelist} ; do
   
         file_yearly=${file/_mon.special.nc/_year.special.nc}   
   
         if [ ! -e ${file_yearly} ] ; then cdo -s yearmean ${file} ${file_yearly}; fi
   
      done

      if [ "${filelist}" == "" ]; then

         cdo setyear,${year} -setrtomiss,-9e99,9e99 ensmin.$((${year}-1)).nc ensmin.${year}.nc
         cdo setyear,${year} -setrtomiss,-9e99,9e99 ensmax.$((${year}-1)).nc ensmax.${year}.nc
	 cdo setyear,${year} -setrtomiss,-9e99,9e99 ensmean.$((${year}-1)).nc ensmean.${year}.nc
	 cdo setyear,${year} -setrtomiss,-9e99,9e99 ensstd.$((${year}-1)).nc ensstd.${year}.nc

      else      

         cdo -s ensmin $(echo ${filelist} | sed -e "s/_mon.special.nc/_year.special.nc/g") ensmin.${year}.nc  
         cdo -s ensmax $(echo ${filelist} | sed -e "s/_mon.special.nc/_year.special.nc/g") ensmax.${year}.nc
         cdo -s ensmean $(echo ${filelist} | sed -e "s/_mon.special.nc/_year.special.nc/g") ensmean.${year}.nc
         cdo -s ensstd $(echo ${filelist} | sed -e "s/_mon.special.nc/_year.special.nc/g") ensstd.${year}.nc

      fi

      cdo -s -f nc -setvar,nruns -setyear,${year} -input,${gridref_timeseries} nfiles.${year}.nc << __EOR__   
      ${nfiles}
__EOR__
   
   done # year
   
   cdo -O mergetime ensmin.????.nc ${ensstatdir}/$(echo ${ensstat_filename} | sed -e "s/ENSSTAT/ensmin/g" -e "s/TS/year/g" )
   cdo -O mergetime ensmax.????.nc ${ensstatdir}/$(echo ${ensstat_filename} | sed -e "s/ENSSTAT/ensmax/g" -e "s/TS/year/g" )
   cdo -O mergetime ensmean.????.nc ${ensstatdir}/$(echo ${ensstat_filename} | sed -e "s/ENSSTAT/ensmean/g" -e "s/TS/year/g" )
   cdo -O mergetime ensstd.????.nc ${ensstatdir}/$(echo ${ensstat_filename} | sed -e "s/ENSSTAT/ensstd/g" -e "s/TS/year/g" )
   
   ln -s ${ensstatdir}/$(echo ${ensstat_filename} | sed -e "s/ENSSTAT/ensmin/g" -e "s/TS/year/g") ensmin.allyears.nc
   ln -s ${ensstatdir}/$(echo ${ensstat_filename} | sed -e "s/ENSSTAT/ensmax/g" -e "s/TS/year/g") ensmax.allyears.nc
   ln -s ${ensstatdir}/$(echo ${ensstat_filename} | sed -e "s/ENSSTAT/ensmean/g" -e "s/TS/year/g") ensmean.allyears.nc
   ln -s ${ensstatdir}/$(echo ${ensstat_filename} | sed -e "s/ENSSTAT/ensstd/g" -e "s/TS/year/g") ensstd.allyears.nc
   
   # lyear=$(echo "1 + $(cdo -s showyear ensmin.allyears.nc | awk '{print $NF}')" | bc)
   # fyear=$(cdo -s showyear ensmin.allyears.nc | awk '{print $1}')
   
   outdir=${dir}/monitoring/${set}/plots/timeseries
   if [ ! -e ${outdir} ] ; then mkdir -p ${outdir}; fi
   
   for var in tsurfgm totprecgm rad0gm ; do
   
      outfile=${outdir}/timeseries_annual_${var}.png
   
      cat << __EOF__ | sed -e "s/^ *//" > plot.py
          import xarray as xr
          import numpy as np
          import matplotlib.pyplot as plt
      
          dsmean=xr.open_dataset('ensmean.allyears.nc')
          dsmin=xr.open_dataset('ensmin.allyears.nc')
          dsmax=xr.open_dataset('ensmax.allyears.nc')
          dsstd=xr.open_dataset('ensstd.allyears.nc')
   
          #-- create figure and axes instances
          dpi = 100
          fig = plt.figure(figsize=(1100/dpi, 1100/dpi), dpi=dpi)
          ax  = fig.add_axes([0.1,0.1,0.8,0.9])
          x=dsmean.time.values
   
          mean=dsmean.${var}.isel(lon=0,lat=0)
          std=dsstd.${var}.isel(lon=0,lat=0)
          spreadmax=mean+std
          spreadmin=mean-std
          plt.plot(x,mean,linewidth=4,color='red')
          # plt.plot(x,spreadmin,linestyle='dashed',linewidth=4,color='red')
          # plt.plot(x,spreadmax,linestyle='dashed',linewidth=4,color='red')
   
          ax.fill_between(x,dsmin.${var}.isel(lon=0,lat=0),dsmax.${var}.isel(lon=0,lat=0),facecolor='gray')
   
          plt.xlabel('year', fontsize=30)
          plt.ylabel('${var}', fontsize=30)
    
          # locs, labels=plt.xticks()
          # print(locs)
    
          ax.tick_params(axis='both', labelsize=20)
   
          plt.savefig('${outfile}', bbox_inches='tight', dpi=dpi)
      
__EOF__
      
      python3 plot.py
   
   done # var

   cd ${oldpwd}
   # rm -rf ${tmp}

done # set

exit 0
