function interpolate_SST() {
  local DYNAMICS_GRID_FILE="$1"
  local SST_INPUT="$2"
  local SST_OUTPUT="$3"

  remap_namelist="tmp_sst_remap.nmp"
  cat >${remap_namelist} <<END

  &remap_nml
  in_filename = "${SST_INPUT}"
  in_type = 1
  out_grid_filename = "${DYNAMICS_GRID_FILE}"
  out_filename = "${SST_OUTPUT}"
  out_type = 2
  out_filetype = 4
  /

  &input_field_nml
  inputname = "T_SEA"
  outputname = "T_SEA"
  code=167
  /

END

  ulimit -s unlimited
  OMP_NUM_THREADS=1 iconremap --remap_nml ${remap_namelist} -vvv

  # Remove namelist
  rm -f ${remap_namelist}

}

#!/bin/bash

# Define the function
function grid_filename_from_extpar() {
  # Assuming the filename is given as the first argument to the function
  filename="$1"

  # This extracts the grid identifier from the filename
  grid_id="${filename:12:13}"

  # This constructs the new filename
  grid_filename="icon_grid_${grid_id}.nc"

  # This prints out the new filename
  echo "${grid_filename}"
}

function interpolate_extpar() {

  local DYNAMICS_GRID_FILE="$1"
  local EXTERNAL_PARAMETERS_FILE="$2"

  echo "Converting extpar file to ICON grid"

  EXTERNAL_PARAMETERS_GRID_FILE=$(grid_filename_from_extpar ${EXTERNAL_PARAMETERS_FILE})

  if [ "${DYNAMICS_GRID_FILE}" != "${EXTERNAL_PARAMETERS_GRID_FILE}" ]; then
    cat >remap.nml <<END
&remap_nml
in_filename = "${EXTERNAL_PARAMETERS_FILE}"
in_grid_filename = "${EXTERNAL_PARAMETERS_GRID_FILE}"
in_type = 2
out_grid_filename = "${DYNAMICS_GRID_FILE}"
out_filename = "tmp_extpar.nc"
out_type = 2
out_filetype = 4
/
END

    cat <<END >extend_remap_namelist.py
import xarray as xr

with xr.open_dataset("${EXTERNAL_PARAMETERS_FILE}") as ds, open('remap.nml', 'a') as file:
   for ida, da in enumerate(ds):
       if da in ["clat","clon","lat","lon"]:
           continue
       file.write('&input_field_nml\n')
       file.write(f'inputname = "{da}"\n')
       file.write(f'outputname = "{da}"\n')
       file.write('/\n\n')
END

    python3 extend_remap_namelist.py

    ulimit -s unlimited
    OMP_NUM_THREADS=1 iconremap --remap_nml remap.nml #-vvv
    rm remap.nml

    # Somehow, when we interpolate we lose some attributes.
    # Here we solve this issue adding the attributes.
    cat <<END >add_attributes.py
# Add sst field to analysis file.
from enstools.io import read

# Open files
with read("tmp_extpar.nc") as ds, read("${EXTERNAL_PARAMETERS_FILE}") as original_extpar_ds:
    if "rawdata" in original_extpar_ds.attrs:
        ds.attrs["rawdata"] = original_extpar_ds.attrs["rawdata"]
        ds.to_netcdf("extpar.nc")
END
    python3 add_attributes.py

    rm "tmp_extpar.nc"

  else
    ln -sf "${EXTERNAL_PARAMETERS_FILE}" "extpar.nc"
  fi
}

function interpolate_analysis() {

  local INPUT_GRID="$1"
  local ANALYSIS_FILE="$2"
  local OUTPUT_GRID="$3"

  echo "Converting analysis to netcdf"

  cat >remap.nml <<END
&remap_nml
in_filename = "${ANALYSIS_FILE}"
in_grid_filename = "${INPUT_GRID}"
in_type = 2
out_grid_filename = "${OUTPUT_GRID}"
out_filename = "interpolated_analysis.nc"
out_type = 2
out_filetype = 4
/
END

  cat <<END >extend_remap_namelist.py
from enstools.io import read

with read("${ANALYSIS_FILE}") as ds, open('remap.nml', 'a') as file:
    all_vars = list(ds.data_vars)
    remap_vars = []
    for var in all_vars:
        if not "bnds" in var and not '_vertices' in var and not 'lat' in var and not 'lon' in var:
            remap_vars.append(var)

    for var in remap_vars:
        file.write('&input_field_nml\n')
        file.write(f'inputname = "{var}"\n')
        file.write(f'outputname = "{var}"\n')
        file.write('/\n\n')
END

  python3 extend_remap_namelist.py

  ulimit -s unlimited
  OMP_NUM_THREADS=1 iconremap --remap_nml remap.nml #-vvv
  rm remap.nml
}

function integrate_sst_to_analysis() {
  local INTERPOLATED_SST="$1"
  local ANALYSIS_FILE="interpolated_analysis.nc"

  cat <<END >integrate_sst_to_analysis.py
# Add sst field to analysis file.
from enstools.io import read
import os
import numpy as np

os.environ['ECCODES_GRIB_NO_INDEX'] = '1'
# Open files
with read("${ANALYSIS_FILE}") as ds, read("${INTERPOLATED_SST}") as sst_ds:
  date=ds.time.values
  day_of_year = int((date - date.astype('datetime64[Y]')) / np.timedelta64(1, 'D'))
   Replace analysis T_SEA with info from the sst_clim file.
  ds["T_SEA"]= sst_ds["T_SEA"].isel(time=day_of_year+1)
  ds.to_netcdf("analysis.nc")
END

  python integrate_sst_to_analysis.py

}

function integrate_sst_to_extpar() {
  local INTERPOLATED_SST="$1"
  local EXTERNAL_PARAMETERS_FILE="$2"

  cat <<END >integrate_sst_to_extpar.py
# Add sst field to analysis file.
from enstools.io import read
import os
import numpy as np
from datetime import datetime, timedelta
import xarray as xr

os.environ['ECCODES_GRIB_NO_INDEX'] = '1'
# Open files
with read("${EXTERNAL_PARAMETERS_FILE}") as ds, read("${INTERPOLATED_SST}") as sst_ds:
  da = sst_ds["T_SEA"]
  year = 2000 # or any other year
  dates = [(datetime(year, 1, 1) + timedelta(days=idx)).month for idx, _ in enumerate(da['time'].values)]

  da['time'] = dates

  # Replace analysis T_SEA with monthly mean from the sst_ds.
  monthly_means = [da.sel(time=month).mean(dim="time") for month in range(1,13)]
  concat_da = xr.concat(monthly_means, dim='time')
  concat_da['time'] = np.arange(1, 13) # replace 'time' with month numbers
  concat_da = concat_da.drop(['clon', 'clat'])

  # Replace analysis T_SEA with monthly mean from the sst_ds.
  ds["T_SEA"] = concat_da

  ds.to_netcdf("extpar.nc")

END

  python integrate_sst_to_extpar.py

}
