#!/bin/bash
set -u

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

# Inputs
in_folder="/storepelican/dverfail/bsc/anom_bias"
oras5_file="ORAS5_preinterp_opa1_grid_T_1985-2014"

# Calculate climatological ORAS5 reference
module --force purge
module load releases/2020b
module load noarch
module load NCO
module load CDO

cd ${in_folder}/"oras5_clim"

cdo ymonmean ${oras5_file}.nc ${oras5_file}_Mm.nc

# Extract January 1985 from ORAS5 file and month of January from climatological reference
ncks -d time_counter,0 ${oras5_file}.nc ORAS5_preinterp_opa1_grid_T_Jan1985.nc
ncks -d time_counter,0  ${oras5_file}_Mm.nc ${oras5_file}_Mm_Jan.nc

# Calculate anomalies: Jan 1985 -climatological ref (these anomalies will then have to be ADDED to the PARASO restarts)
ncdiff -v conservative_temperature ORAS5_preinterp_opa1_grid_T_Jan1985.nc ${oras5_file}_Mm_Jan.nc ../anomalies/thetao_anomaly_3D_Jan_1985.nc
ncdiff -v absolute_salinity ORAS5_preinterp_opa1_grid_T_Jan1985.nc ${oras5_file}_Mm_Jan.nc ../anomalies/so_anomaly_3D_Jan_1985.nc

# Extract anomalies for the surface only
cd ${in_folder}/anomalies
ncks -d deptht,0 thetao_anomaly_3D_Jan_1985.nc tos_anomaly_2D_Jan_1985.nc
ncks -d deptht,0 so_anomaly_3D_Jan_1985.nc sos_anomaly_2D_Jan_1985.nc
ncwa -O -a deptht tos_anomaly_2D_Jan_1985.nc tos_anomaly_2D_Jan_1985.nc
ncwa -O -a deptht sos_anomaly_2D_Jan_1985.nc sos_anomaly_2D_Jan_1985.nc
ncks -O -x -v deptht tos_anomaly_2D_Jan_1985.nc tos_anomaly_2D_Jan_1985.nc
ncks -O -x -v deptht sos_anomaly_2D_Jan_1985.nc sos_anomaly_2D_Jan_1985.nc

# Reformat anomaly files to have the same format as the PARASO restart files
ncrename -O -d time_counter,time -v lat,nav_lat -v lon,nav_lon -v conservative_temperature,sst_m tos_anomaly_2D_Jan_1985.nc
ncrename -O -d time_counter,time -v lat,nav_lat -v lon,nav_lon -v absolute_salinity,sss_m sos_anomaly_2D_Jan_1985.nc
ncrename -O -d time_counter,time -d deptht,z -v deptht,nav_lev -v lat,nav_lat -v lon,nav_lon -v conservative_temperature,tb thetao_anomaly_3D_Jan_1985.nc
ncrename -O -d time_counter,time -d deptht,z -v deptht,nav_lev -v lat,nav_lat -v lon,nav_lon -v absolute_salinity,sb so_anomaly_3D_Jan_1985.nc
cp tos_anomaly_2D_Jan_1985.nc tos_anomaly_2D_Jan_1985_reformatted.nc
cp sos_anomaly_2D_Jan_1985.nc sos_anomaly_2D_Jan_1985_reformatted.nc
ncap2 -s 'tn=tb' thetao_anomaly_3D_Jan_1985.nc thetao_anomaly_3D_Jan_1985_reformatted.nc
ncap2 -s 'sn=sb' so_anomaly_3D_Jan_1985.nc so_anomaly_3D_Jan_1985_reformatted.nc

# Combine all anomaly files into a single one
cp so_anomaly_3D_Jan_1985_reformatted.nc all_anomalies_Jan_1985.nc
ncks -A tos_anomaly_2D_Jan_1985_reformatted.nc sos_anomaly_2D_Jan_1985_reformatted.nc
ncks -A sos_anomaly_2D_Jan_1985_reformatted.nc thetao_anomaly_3D_Jan_1985_reformatted.nc
ncks -A thetao_anomaly_3D_Jan_1985_reformatted.nc all_anomalies_Jan_1985.nc
ncap2 -s 'nav_lat=float(nav_lat); lat_bnds=float(lat_bnds); nav_lon=float(nav_lon); lon_bnds=float(lon_bnds); sn=double(sn); sb=double(sb); sss_m=double(sss_m); sst_m=double(sst_m); tb=double(tb); tn=double(tn);' all_anomalies_Jan_1985.nc all_anomalies_Jan_1985_reformatted.nc