= Interpolation with CDO

[[PageOutline(2)]]

'''You can find here an example of bash script to interpolate forcing core2 files.[[BR]]

To interpolate files you need to do 3 steps:[[BR]]

1. have good grid description file (if not you can create it, see instructions 1.) below)[[BR]]

2. prepare files to be interpolated ((see 2.) below)[[BR]]

3. remapping on ORCA2 grid ((see 3.) below)'''[[BR]]


NOTE: This script works with cdo version 1.4.0.1



== 1) Prepare grid description file for the ORCA curvilinear grid

You need to prepare a netCDF grid description file for the ORCA curvilinear grid.[[BR]]
This file must contain:[[BR]]
    - a longitude variable with exactly 2 dimensions and the attribute: units = "degrees_east"[[BR]]
    - a latitude variable with exactly 2 dimensions and the attribute: units = "degrees_north"[[BR]]
    - a dummy variable with the attribute: coordinates = "xxx yyy" where xxx and yyy are the name of the longitude and latitude variables mentioned just above[[BR]]

There is an example of a [http://www.unidata.ucar.edu/software/netcdf/ netCDF] grid description file for the T grid of ORCA2:[[BR]]


{{{
netcdf grid_ORCA2_T {
dimensions:
        y = 149 ;
        x = 182 ;
variables:
        byte dummy(y, x) ;
                dummy:coordinates = "glamt gphit" ;
        double glamt(y, x) ;
                glamt:units = "degrees_east" ;
        double gphit(y, x) ;
                gphit:units = "degrees_north" ;
}
}}}

There is an example of nco commands to create this file from coordinates.nc (input file used by NEMO for its grid description)[[BR]]


{{{#!sh
# select glamt and gphit variables from the coordinates file
ncks -O -C -a -v glamt,gphit coordinates_ORCA_R2.nc grid_ORCA2_T.nc

}}}

make sure that coordinates variables contains only 2 dimension. Use ncwa -a to remove degenerated dimensions (with a size of 1) for example:[[BR]]


{{{#!sh
# remove degenerated dimention time (if existing)
ncwa -O -a time grid_ORCA2_T.nc  grid_ORCA2_T.nc
# remove degenerated dimention z (if existing)
ncwa -O -a z grid_ORCA2_T.nc  grid_ORCA2_T.nc
}}}

{{{#!sh
# add a dummy variable
ncap2 -O -s 'dummy[y,x]=1b' grid_ORCA2_T.nc grid_ORCA2_T.nc
# add needed attributes
ncatted -a coordinates,dummy,c,c,'glamt gphit' \
        -a units,glamt,c,c,'degrees_east'      \
        -a units,gphit,c,c,'degrees_north'  grid_ORCA2_T.nc
}}}


== 2) Prepare file to be interpolated

for a curvilinear grid:[[BR]]

    - you have to follow above spefication[[BR]]

For a regular lon/lat grid:[[BR]]

    - The coordinate variables must have the same name as the dimensions. You can chose any name. If you use the name x for the dimension then you have to use also x for the coordinate variable.[[BR]]
    - The coordinate variable must also have the attribute units = "degrees_east" or "degrees_north" to identify the lon/lat variables:[[BR]]

{{{
netcdf etopo5 {
 dimensions:
         x = 4320 ;
         y = 2161 ;
 variables:
         float x(x) ;
                  x:units = "degree_east" ;
         float y(y) ;
                 y:units = "degree_north" ;
         float bath(y, x) ;
 }
}}}


== 3) Remapping on ORCA2

'''Generate weights with bilinear interpolation, and then remapping for a scalar field:'''[[BR]]

  
{{{#!sh

    cdo genbil,grid_ORCA2_T.nc $file bil_orca2_weights.nc
    cdo remap,grid_ORCA2_T.nc,bil_orca2_weights.nc $file_fill.nc $file_orca2.nc
}}}


where[[BR]]

  - $file_fill.nc is a file whose missing values is filled (see below interpolation_loop.sh)[[BR]]

  - $file_orca2.nc is the output[[BR]]

  - and bil_orca2_weights.nc is a weight file generated (see below)[[BR]]


''Generate weights with bilinear interpolation, and then remapping for a vectorial field:'''[[BR]]


{{{#!sh

    cdo genbic,grid_ORCA2_T.nc $file bic_orca2_weights.nc
    cdo remap,grid_ORCA2_T.nc,bic_orca2_weights.nc $file_fill.nc $file_orca2.nc
}}}


where[[BR]]

  - $file_fill.nc is a file whose missing values is filled (see below interpolation_loop.sh)[[BR]]

  - $file_orca2.nc is the output[[BR]]

  - and bic_orca2_weights.nc is a weight file generated (see below)[[BR]]


== Example of "interpolation_loop.sh" script

This script needs to have input directory in your execution's directory to copy *.nc files and it calls 4 others sub-scripts: "create_orca2.sh", "prepare_mask.sh", "bilweights.sh" and "bicweights.sh"[[BR]]


NOTE: Before running "interpolation_loop.sh", you nedd to checks different things:[[BR]]

1. control latitude of your mask file. If it goes from north to south and latitude in your file goes from south to north you have to invert it:[[BR]]

{{{
cdo invertlat ${file_mask_land}.nc ${file_mask_land}.SN.nc
}}}


2. control format of land variable, if it is short and your file has float format you will have problems with cdo created missing values, and filling missing values

{{{#!sh
#convert in float (from short) land variable
ncap -O -s "land=float(land)" ${file_mask_land}.SN.nc ${file_mask_land}.SN.float.nc
}}}

3. control values of your land variable, if it is 1 on earth and 0 on ocean is not good if you are interpolating files for ocean forcing

{{{#!sh
#convert land=0 on earth, 1 on ocean
ncap -O -s 'land=1-land' mask.nc mask.nc
}}}

{{{#!sh
#rename "land" in "mask"
ncrename -v land,mask mask.nc
}}}


4. control values of calendar's attribute in your file, if it is "NOLEAP" it is better to convert it in a recognized format by nccated (ex. 365_day)

{{{#!sh
ncatted -a calendar,${variable},m,c,365_day ${file}.nc ${file}_365.nc
}}}

 
'''main programm: interpolation_loop.sh '''

{{{#!sh
#!/bin/bash
#script to interpolate files on ORCA2 grids

set -u
LC_ALL=C; export LC_ALL

type cdo
status_cdo=${?}
if [ ${status_cdo} != 0 ] ; then
        echo "ERROR: cdo not found in your PATH"
        exit 1
fi

file_orca2="grid_ORCA2_T.nc"
#build ORCA2_grid if necessary
if [ ! -f ${file_orca2} ] ; then
 /bin/sh create_orca2.sh
fi

#preparing mask file if necessary
if [ ! -f mask.nc ] ; then
 /bin/sh prepare_mask.sh
fi

dir=/PATH/FILES/TO/BE/INTERPOLATED

`cp ${dir}/*.nc input/.`
cd input

for file in $(echo `ls *.nc | sed -e "s/\.nc//"`) ; do
 
 # create missing values over land with e.g. 'ifthen':
 echo "create missing values ifthen mask"
 cdo ifthen ../mask.nc ${file}_365.nc ${file}_mask.nc
 status_ifthen=${?}
 if [ ${status_ifthen} != 0 ] ; then
         echo "ERROR ifthen mask"
         exit 1
 fi
 
 
 # fill missing value (here over land):
 echo "fill missing values"
 cdo fillmiss ${file}_mask.nc ${file}_fill.nc
 status_fillmiss=${?}
 if [ ${status_fillmiss} != 0 ] ; then
  echo "ERRORE cdo fillmiss"
         exit 1
 fi
 
 
 #remapping on ORCA2 grid
 echo "remapping "
 cd ..
# looking for first part of filename, to see if it is vectorial field, for interannual files
 name=$(echo $file | awk -F \_ {'print $1'} )
 if [ ${name} = "u" -o ${name} = "v" ] ; then
    echo "vectorial field"
    /bin/sh bicweights.sh input/${file}_fill.nc
    cdo remap,grid_ORCA2_T.nc,bic_orca2_weights.nc input/${file}_fill.nc input/${file}_orca2.nc
 else 
    echo "scalar field"
    /bin/sh bilweights.sh input/${file}_fill.nc
    cdo remap,grid_ORCA2_T.nc,bil_orca2_weights.nc input/${file}_fill.nc input/${file}_orca2.nc
 fi
 status_remap=${?}
 if [ ${status_remap} != 0 ] ; then
         echo "ERROR cdo remap"
         exit 1
 fi
 echo "status remap: ${status_remap}"
 cd input
done
}}}

'''first sub-programm: file create_orca2.sh'''
{{{#!sh
#!/bin/bash
#script to create ORCA2 grid

set -u
LC_ALL=C ; export LC_ALL

echo "create orca2 grid"
ncks -O -C -a -v glamt,gphit coordinates.nc grid_ORCA2_T.nc 

#remove not necessar axes:
ncwa -O -a time,z grid_ORCA2_T.nc grid_ORCA2_T.nc

#add dummy variable:
ncap -O -s 'dummy[y,x]=1b' grid_ORCA2_T.nc grid_ORCA2_T.nc

#add needed attributes
ncatted -a coordinates,dummy,c,c,'glamt gphit'  \
        -a units,glamt,c,c,'degrees_east'       \
        -a units,gphit,c,c,'degrees_north' grid_ORCA2_T.nc
}}}

'''Second sub-programm: prepare_mask.sh'''
{{{#!sh
#!/bin/bash
#preparing mask file

set -u
LC_ALL=C ; export LC_ALL

file_mask_land="land.sfc.gauss"
# invert latitude (from SOUTH to NORTH)
cdo invertlat ${file_mask_land}.nc ${file_mask_land}.SN.nc

#convert in float (from short) land variable
ncap -O -s "land=float(land)" ${file_mask_land}.SN.nc ${file_mask_land}.SN.float.nc 

# remove all attribute to land
ncatted -a ,land,d,, ${file_mask_land}.SN.float.nc ${file_mask_land}.SN.float.miss.nc

# extract land
#ncks -O -v land land.sfc.gauss.SN.float.nc mask.nc 
ncks -O -v land ${file_mask_land}.SN.float.miss.nc mask.nc 
ncwa -O -a time mask.nc mask.nc

#convert land=0 on earth, 1 on ocean
ncap -O -s 'land=1-land' mask.nc mask.nc

#rename "land" in "mask"
ncrename -v land,mask mask.nc
chmod 755 mask.nc
}}}


'''Third sub-programm: bicweights.sh'''
{{{#!sh
#!/bin/bash

#remapping on ORCA2 grid
echo "gen bicubic WEIGHTS"
export REMAP_EXTRAPOLATE=on
cdo genbic,grid_ORCA2_T.nc $1 bic_orca2_weights.nc
}}}


'''Forth sub-programm: bicweights.sh'''
{{{#!sh
#!/bin/bash

#remapping on ORCA2 grid
echo "gen bilinear WEIGHTS"
export REMAP_EXTRAPOLATE=on
cdo genbil,grid_ORCA2_T.nc $1 bil_orca2_weights.nc
}}}