#!/bin/bash
set -u

# This script takes EC-Earth hindcast outputs, reinterpolates them to the NEMO grid, and then generates BDY dta from it.
# Author: C. Pelletier & D. Verfaillie

# Inclusive year range
year_start=1985
year_stop=2014

# Input prefix (with folder) for input data files.
in_prefix="/storepelican/dverfail/bsc/anom_bias/hindcasts"
# Name of experience
name_exp="a3du"
# Frequency
freq="1m"

#for fc_num in fc00 fc01 fc02; do
for fc_num in fc00 fc01 fc02; do

# Working folder where pre-interpolated datafiles will be put
dest_folder="/storepelican/dverfail/bsc/nemo_bdy/bdydta"
dest_prefix="ECEarth_${name_exp}_${fc_num}_${freq}_eORCA025-SO_BDYDTA_nbr1_"

pwd=`pwd`
# Python script for reformatting BDYDTA (final step)
reformat_script="/home/elic/dverfail/PARAMOUR_TOOLS/nemo_bdy/generate_new_bdy_dta_onefile.py"

# min_lon, max_lon, min_lat, max_lat
# To precut the input data prior to interpolation (not required, but less expensive if your source code is much bigger than the destination one, which is the case from global to the Southern Ocean.
# Take a generous margin (4-5 nodes bigger than the destination grid limit)
latlon_cut="-180,180,-90,-28"

# # Script for converting EOS80 to TEOS-10 (look at the script comments for why it exists and what it does).
# not needed for EC-Earth, which is already in TEOS-10
# curr_dir=`pwd`
# eos10_script="${curr_dir}/convert_eos10_new.py"

# NEMO destination CDO grid description (see `cdo_get_griddes.sh`).
dest_gridfile_t="/storepelican/pelletie/nemo_bdy/grids/eORCA025-SO_grid_T.grid"
dest_gridfile_u="/storepelican/pelletie/nemo_bdy/grids/eORCA025-SO_grid_U.grid"
dest_gridfile_v="/storepelican/pelletie/nemo_bdy/grids/eORCA025-SO_grid_V.grid"

pid="$$"

mkdir -p ${dest_folder}
cd ${dest_folder}

for((year=year_start;year<=year_stop;year++)); do

    # Loop over input years.
    
    module --force purge
    module load releases/2020b
    module load noarch
    module load NCO
    module load CDO

    # 1) Extract T, S and ssh (T grid).

    cdo selname,thetao,so ${in_prefix}/${fc_num}/${name_exp}_${freq}_opa_grid_T_3D_${year}.nc tmp_pid${pid}.nc
    cdo selname,zos ${in_prefix}/${fc_num}/${name_exp}_${freq}_opa_grid_T_2D_${year}.nc tmp2_pid${pid}.nc
    cdo -O merge tmp_pid${pid}.nc tmp2_pid${pid}.nc tmp3_pid${pid}.nc

    # precut
    cdo -O sellonlatbox,${latlon_cut} tmp3_pid${pid}.nc tmp_pid${pid}.nc

    ncrename -v olevel,deptht tmp_pid${pid}.nc
    ncrename -d olevel,deptht tmp_pid${pid}.nc
    
    # remap
    cdo remapbil,${dest_gridfile_t} tmp_pid${pid}.nc tmp2_pid${pid}.nc
    cdo setmisstonn tmp2_pid${pid}.nc tmp_pid${pid}.nc

    module --force purge
    module load releases/2019a
    module load noarch
    module load Java/11.0.2
    module load ELIC_Python/1-foss-2019a-Python-3.7.2

    input_full=`readlink -f tmp_pid${pid}.nc`
    output_full="${dest_folder}/${dest_prefix}grid_T_y${year}.nc"
    
    python ${reformat_script} ${input_full} t ${output_full}

    module --force purge
    module load releases/2020b
    module load noarch
    module load NCO
    module load CDO

    # 2) Extract u velocities.

    #calculate ubarotropic from uo (3D): vertical mean of uo
    cdo selname,uo ${in_prefix}/${fc_num}/${name_exp}_${freq}_opa_grid_U_3D_${year}.nc tmp_pid${pid}.nc
    cdo vertmean tmp_pid${pid}.nc tmp2_pid${pid}.nc
    ncap2 -O -s 'ubar = uo;' tmp2_pid${pid}.nc tmp3_pid${pid}.nc
    cdo -O merge tmp_pid${pid}.nc tmp3_pid${pid}.nc tmp4_pid${pid}.nc

    # precut
    cdo -O sellonlatbox,${latlon_cut} tmp4_pid${pid}.nc tmp_pid${pid}.nc

    ncap2 -O -s 'uobaroclinic = uo - ubar;' tmp_pid${pid}.nc tmp2_pid${pid}.nc
    ncks -x -v uo,uo_2 tmp2_pid${pid}.nc -O tmp_pid${pid}.nc
    ncrename -v ubar,uobarotropic tmp_pid${pid}.nc
    ncrename -v olevel,depthu -v olevel_bnds,depthu_bnds tmp_pid${pid}.nc
    ncatted -a bounds,depthu,o,c,"depthu_bnds" -a name,depthu,o,c,"depthu" tmp_pid${pid}.nc
    ncrename -d olevel,depthu tmp_pid${pid}.nc
    ncpdq -O -a time_counter,depthu,y,x tmp_pid${pid}.nc tmp2_pid${pid}.nc
    
    # remap
    cdo remapbil,${dest_gridfile_u} tmp2_pid${pid}.nc tmp_pid${pid}.nc
    cdo setmisstonn tmp_pid${pid}.nc tmp2_pid${pid}.nc

    module --force purge
    module load releases/2019a
    module load noarch
    module load Java/11.0.2
    module load ELIC_Python/1-foss-2019a-Python-3.7.2

    input_full=`readlink -f tmp2_pid${pid}.nc`
    output_full="${dest_folder}/${dest_prefix}grid_U_y${year}.nc"
    
    python ${reformat_script} ${input_full} u ${output_full}

    module --force purge
    module load releases/2020b
    module load noarch
    module load NCO
    module load CDO
    
    # 3) Extract v velocities.

    #calculate vbarotropic from vo (3D): vertical mean of vo
    cdo selname,vo ${in_prefix}/${fc_num}/${name_exp}_${freq}_opa_grid_V_3D_${year}.nc tmp_pid${pid}.nc
    cdo vertmean tmp_pid${pid}.nc tmp2_pid${pid}.nc
    ncap2 -O -s 'vbar = vo;' tmp2_pid${pid}.nc tmp3_pid${pid}.nc
    cdo -O merge tmp_pid${pid}.nc tmp3_pid${pid}.nc tmp4_pid${pid}.nc

    # precut
    cdo -O sellonlatbox,${latlon_cut} tmp4_pid${pid}.nc tmp_pid${pid}.nc

    ncap2 -O -s 'vobaroclinic = vo - vbar;' tmp_pid${pid}.nc tmp2_pid${pid}.nc
    ncks -x -v vo,vo_2 tmp2_pid${pid}.nc -O tmp_pid${pid}.nc
    ncrename -v vbar,vobarotropic tmp_pid${pid}.nc
    ncrename -v olevel,depthv -v olevel_bnds,depthv_bnds tmp_pid${pid}.nc
    ncatted -a bounds,depthv,o,c,"depthv_bnds" -a name,depthv,o,c,"depthv" tmp_pid${pid}.nc
    ncrename -d olevel,depthv tmp_pid${pid}.nc
    ncpdq -O -a time_counter,depthv,y,x tmp_pid${pid}.nc tmp2_pid${pid}.nc
    
    # remap
    cdo remapbil,${dest_gridfile_v} tmp2_pid${pid}.nc tmp_pid${pid}.nc
    cdo setmisstonn tmp_pid${pid}.nc tmp2_pid${pid}.nc 

    module --force purge
    module load releases/2019a
    module load noarch
    module load Java/11.0.2
    module load ELIC_Python/1-foss-2019a-Python-3.7.2

    input_full=`readlink -f tmp2_pid${pid}.nc`
    output_full="${dest_folder}/${dest_prefix}grid_V_y${year}.nc"
    
    python ${reformat_script} ${input_full} v ${output_full}

    module --force purge
    module load releases/2020b
    module load noarch
    module load NCO
    module load CDO
    
done
done

rm -rf ./*pid${pid}*