#!/bin/bash
set -u

# This script generates CDO grid file descriptions.
# These CDO description files are used for horizontal interpolation from/to a given grid.
# The script assumes that the horizontal dimensions are called (y,x) but feel free to edit.

# Basically, from the netCDF gridfile, this scripts creates a dummy netCDF that CDO recognizes, to get the CDO grid description file. We'll do that:
# - longitude variable with netCDF attributes standard_name="longitude", units="degrees East"
# - latitude variable with netCDF attributes standard_name="latitude", units="degrees North"
# - one dummy variable with a netCDF attribute coordinates='<lat> <lon>'

# Then cdo griddes dummy.nc > my_grid.grid will work.

# Input file with lat, lon and depth
input_file="/mfast/pelletie/oras5/grids/mesh_mask.nc"
lat_var="gphit"
lon_var="glamt"
depth_var="gdept_1d"

output_griddes="/mfast/pelletie/oras5/grids/example_mine.grid"

module purge
module load 2018b CDO NCO

pid="$$"

# Just take lat/lon from file
ncks -C -v ${lat_var},${lon_var} ${input_file} -O tmp_pid${pid}.nc
# Edit lat attributes
ncatted -a standard_name,${lat_var},o,c,"latitude" -a units,${lat_var},o,c,"degrees North" tmp_pid${pid}.nc
# Edit lon attributes
ncatted -a standard_name,${lon_var},o,c,"longitude" -a units,${lon_var},o,c,"degrees East" tmp_pid${pid}.nc
# Create some dummy variable. Here, dimension names are assumed to be y/x. Change if needed.
ncap2 -O -s "dummy[\$y,\$x]=1.f;" tmp_pid${pid}.nc tmp2_pid${pid}.nc
# Edit dummy variable attributes to say its coordinates are you lat/lons.
ncatted -a coordinates,dummy,o,c,"${lat_var} ${lon_var}" tmp2_pid${pid}.nc

# Get CDO griddes
cdo griddes tmp2_pid${pid}.nc > ${output_griddes}

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



# Expert level: if you want to perform conservative interpolation, you'll need grid corners as well.
# For ORCA grids (the NEMO grids), this can be done using F nodes (grid corners), check `get_orca_corners.py`
# Then, you need:
# 1) One additional dimension <corn> with size 4 (the name does not matter) and two additional fields:
#   1a) <lon_corn>(y,x,<corn>): longitudes of grid corners
#   1b)  <lat_corn>(y,x,<corn>): latitudes of grid corners
# 2) Add an attributes "bounds" to the lat/lon fields in the dummy file, i.e.:
#     ncatted -a standard_name,${lat_var},o,c,"latitude" -a units,${lat_var},o,c,"degrees North" -a bounds,${lat_var},o,c,<lat_corn> tmp_pid${pid}.nc
#    ncatted -a standard_name,${lon_var},o,c,"longitude" -a units,${lon_var},o,c,"degrees East" -a bounds,${lon_var},o,c,<lon_corn> tmp_pid${pid}.nc
# Then, create the dummy var the same way as above. The CDO gridfile description should then be much larger (containing corners). And can be used for conservative interpolation.