#!/bin/bash
#
# -- Author : François Massonnet, francois.massonnet@ic3.cat
# -- Date   : 30 Jan 2015
# -- At     : IC3, Barcelona
# -- Modified : 19 Jan 2016, omar.bellprat@bsc.es 
#
# -- Purpose: Generation of an arbitrary number of NEMO oceanic restarts that are copies of a reference, plus a perturbation
#
# -- Method : The reference file is duplicated in this script, then read by a R script. The perturbation is introduced and finally written to this latter file.
#             The script must be placed into the restart directory (e.g. NEMO_Restart_23). The generated restarts have to be renamed after generation. This 
#             script has been tested on MareNostrum3.
#
# -- Input  : NEMO ocean restart from an EC-Earth 3.1 run
# -- Output : N restarts with the same name,but with an index fc0, fc1, ... fcN-1 appended
#
# -- Limitations: Only the surface conditions are perturbed (level index: 1) but this can be changed in the R script

set -o errexit
set -o nounset
set -x

if [ $# == 0 ] ; then
  echo "gener_perturbation.bash ocean_restart_file Nmembers"
  exit
fi

filein=$1
nmemb=$2

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

var=tn          # Variable to be perturbed
per=0.0001      # Standard deviation of gaussian perturbation to be applied,
                # in units of the variable (for tn: in K for example)

for jmemb in `seq 0 $(( $nmemb -1 ))`
do
  echo $jmemb
  # 1. Make a copy of the original file, with the new name
  filenew="${filein%.nc}_fc${jmemb}.nc"
  cp $filein ${filein}.backup
  cp $filein $filenew

  # 2. Prepare the R script

  echo "#!/usr/bin/env Rscript
  library(ncdf)

  # François Massonnet, 30 Jan 2015
  # Adds a gaussian perturbation at the first level of a 3D field
  # Tested only for NEMO restarts
  #
  # This script should be called by a bash script so that the variable and file names are specified, as well as the perturbation

  varname='$var'
  filein <- '$filenew'
  ex.nc           <- open.ncdf(filein,write=TRUE)
  spert <- $per

  myvar     <- get.var.ncdf(ex.nc, varname)
  myvarpert <- myvar
  for (i in seq(1,dim(myvar)[1])){
    for (j in seq(1,dim(myvar)[2])){
      if (myvar[i,j,1] != 0){
        myvarpert[i,j,1] = myvarpert[i,j,1] + rnorm(1,sd=spert)
      }
    }
  }

  put.var.ncdf(ex.nc,varname,myvarpert)
  close.ncdf(ex.nc)" > Rtmpfile.R

  chmod 744 Rtmpfile.R

  # 3. Run the R script, that produces the new NetCDF
  ./Rtmpfile.R 
done