#!/bin/bash
set -u

# This script takes ORAS5 and EC-Earth historical ocean temperature and salinity conditions and calculates anomalies from them.
# Author: D. Verfaillie

# Inputs
in_folder="/storepelican/dverfail/bsc/anom"
oras5_file="ORAS5_preinterp_opa1_grid_T_1985-2014"
histo_1m_thetao="thetao_Omon_EC-Earth3_historical_ensmean_gn_198501-201412"
histo_1m_so="so_Omon_EC-Earth3_historical_ensmean_gn_198501-201412"

# Prepare ORAS5 file to have the same grid and format as the EC-Earth files
module --force purge
module load releases/2020b
module load noarch
module load NCO
module load CDO

cd ${in_folder}/"oras5_clim"

#cdo remapbil,${in_folder}/oras5_clim/ECEarth_grid_T_new.grid ${oras5_file}.nc ${oras5_file}_remapped.nc
#ncrename -O -d time_counter,time -d x,i -d y,j -d nvertex,vertices -d deptht,lev -v time_counter,time -v nav_lon,longitude -v nav_lat,latitude -v nav_lon_bnds,longitude_bnds -v nav_lat_bnds,latitude_bnds -v deptht,lev -v conservative_temperature,thetao -v absolute_salinity,so ${oras5_file}_remapped.nc
#ncap2 -O -s 'defdim("bnds",2);make_bounds(lev,$bnds,"lev_bnds");' ${oras5_file}_remapped.nc ${oras5_file}_remapped_bounds.nc
#ncks -A -v i,j,lev_bnds ../histo_clim/${histo_1m_thetao}.nc ${oras5_file}_remapped_bounds.nc
ncap2 -s 'lev=double(lev); lev_bnds=double(lev_bnds); longitude=float(longitude); longitude_bnds=float(longitude_bnds); latitude=float(latitude); latitude_bnds=float(latitude_bnds);' ${oras5_file}_remapped_bounds.nc ${oras5_file}_remapped_bounds_reformatted.nc
ncks -O -x -v so,sossheig ${oras5_file}_remapped_bounds_reformatted.nc thetao_ORAS5_preinterp_opa1_grid_T_1985-2014.nc
ncks -O -x -v thetao,sossheig ${oras5_file}_remapped_bounds_reformatted.nc so_ORAS5_preinterp_opa1_grid_T_1985-2014.nc

# Calculate monthly and climatological anomalies: histo - oras5 (these anomalies will then have to be SUBSTRACTED from the hindcasts)
ncdiff ../histo_clim/${histo_1m_so}.nc thetao_ORAS5_preinterp_opa1_grid_T_1985-2014.nc ../anom_histo-oras5/thetao_anomaly_1985-2014_1m.nc
ncdiff ../histo_clim/${histo_1m_thetao}.nc so_ORAS5_preinterp_opa1_grid_T_1985-2014.nc ../anom_histo-oras5/so_anomaly_1985-2014_1m.nc

cd ${in_folder}/anom_histo-oras5
cdo ymonmean thetao_anomaly_1985-2014_1m.nc thetao_anomaly_1985-2014_Mm.nc
cdo ymonmean so_anomaly_1985-2014_1m.nc so_anomaly_1985-2014_Mm.nc

# Reformat the anomaly files to have same format as EC-Earth hindcast ocean temperature and salinity files
ncrename -O -d time,time_counter -d i,x -d j,y -d vertices,nvertex -d lev,olevel -d bnds,axis_nbounds -v time,time_counter -v time_bnds,time_counter_bounds -v latitude,nav_lat -v latitude_bnds,bounds_nav_lat -v longitude,nav_lon -v longitude_bnds,bounds_nav_lon -v lev,olevel -v lev_bnds,olevel_bounds thetao_anomaly_1985-2014_Mm.nc thetao_anomaly_1985-2014_Mm_reformatted.nc
ncrename -O -d time,time_counter -d i,x -d j,y -d vertices,nvertex -d lev,olevel -d bnds,axis_nbounds -v time,time_counter -v time_bnds,time_counter_bounds -v latitude,nav_lat -v latitude_bnds,bounds_nav_lat -v longitude,nav_lon -v longitude_bnds,bounds_nav_lon -v lev,olevel -v lev_bnds,olevel_bounds so_anomaly_1985-2014_Mm.nc so_anomaly_1985-2014_Mm_reformatted.nc