SOSIE is Only a Surface Interpolation Environment
About SOSIE

SOSIE is a tool that allows fast and high-quality 2D and 3D interpolation of geophysical fields from a gridded domain to another. "Sosie" is the french word for "double" (like in look-alike). SOSIE is developed and distributed under the GNU General Public License (GPL).
SOSIE was originally intended to interpolate geophysical fields onto the ORCA family of tripolar grids used by the NEMO OGCM (ORCA2, ORCA1, ORCA05, ORCA025). It now supports a wide range of source/target grid configuration for scalar fields interpolation but vector rotation in distorted target grid regions is only supported for the ORCA grids so far.

The main interpolation algorithm of SOSIE is based on the method of Akima (1970): "A New Method of Interpolation and Smooth Surface Fitting Based On Local Procedures, J.of Applied Comput. Math., 17, 589-602." In SOSIE, this method has been coded from scratch in plain Fortran-90. The algorithm is made highly efficient by skipping solving a 16x16 linear system for each treated point of the target grid. Instead, a particularly efficient way to calculate the general solution of the 16x16 system was found (J.M. Brankart, personal communication).

Compared to more widely used interpolation methods such as bilinear or bicubic splines, the Akima method allows, at an extremely low numerical cost, continuous and smooth interpolated fields without errors related to overshoots (as for polynomial functions, see 1D illustrations below).
SOSIE can perform 3D interpolation. This is fake 3D interpolation though: each level is independently interpolated on the target horizontal domain and vertical interpolation (using Akima 1D algorithm) is then performed. SOSIE uses Netcdf as IO file format.

Author: Laurent Brodeau (brodeau (at) gmail)
Contributors: Jean-Michel Brankart, Jean-Marc Molines and Pierre Mathiot

Limitation of the Akima method

The Akima algorithm can only be used for interpolating fields given on non-distorted horizontal grids, i.e. on so-called lat-lon regular domains in which latitude and longitude arrays are 1D and only dependent of respectively j and i. A bilinear interpolation alternative is included in SOSIE and can be used for distorted input domains. However, there is no limitation regarding the type of the target grid: both regular and distorted target grids are supported by the Akima method.

Illustration 1D

Fig. 1-2: Comparison between the Akima method coded in SOSIE and the widely-used cubic spline method (from Matlab), for two different "discontinuous" patterns of input (x,y) points.

Illustration 2D

In this first example, the mean long term sea surface temperature of Reynolds (2002) is interpolated from a spherical input lat-lon grid (360x180) to the tripolar ORCA1 (distorted) grid. This example is provided as a template configuration (example #1) in the SOSIE package, see the "A few examples" section.

Fig. 3: March mean SST (Reynolds, 2002), original regular lat-lon grid (360x180).


Fig. 4: March mean SST (Reynolds, 2002), interpolated on the tripolar ORCA1 grid using the Akima method


In the following example, a given SST snapshot from a random NEMO-ORCA1 simulation is interpolated with the SOSIE bi-linear algorithm from the tripolar ORCA1 (distorted) grid to a spherical output lat-lon grid auto-generated by SOSIE (1x1 degree). This example is provided as a template configuration (example #3) in the SOSIE package, see the "A few examples" section. This example also illustrates the action of the DROWN algorithm, used to extrapolate sea values over continents (see Fig. 6)
Fig. 5: Given NEMO-ORCA1 SST snapshot, on its original ORCA1 grid. Note that not all iso-longitudes and -latitudes are represented for clarity.

Fig. 6: Given NEMO-ORCA1 SST snapshot, interpolated on a regular 1x1 degree lat-lon grid with the bilinear method

Obtaining SOSIE

Up to date version of SOSIE can be installed via Subversion using the following command:

  svn co svn://svn.code.sf.net/p/sosie/code/trunk sosie
        
You can browse the svn depository or download a tarball containing the last svn version of the code on the following page: http://sosie.svn.sourceforge.net/viewvc/sosie/

Using SOSIE

Once you have read the README file and compiled the executable, you simply need to configure everything into a Fortran namelist. Template for the namelist:

!! ------------------- 
!! Namelist for SOSIE 
!! ------------------- 
!! 
!! 
!! *********************** 
!! Input characteristics : 
!! *********************** 
!! 
!! ivect : vector correction control on treated field [integer] 
!!         ivect = 0 : input field is not a component of a vector 
!!         or the target grid is regular (lregout = T) 
!!         * if non-regular distorted target grids (like ORCAX): 
!!         ivect = 1 : input field is a zonal (X) component of a vector 
!!         ivect = 2 : input field is a meridional (Y) component of a vector 
!! 
!! lregin : is the source grid regular? [logical] 
!!          (ie : are input longitude and latitude 1D?) 
!! 
!! cf_in   : file containing the input field to be interpolated [char]
!! cv_in   : name of treated variable (in input field file) [char]
!!
!! cv_t_in : name of time record variable in the input file [char]
!!           or 'missing_rec' if no time record is present in the input file
!!
!! jt1     : first time record to be interpolated
!! jt2     : last  time record to be interpolated
!!           => set jt1 and jt2 to 0 if you want to skip this option 
!!              and interpolate the nt time steps of the current field
!! 
!! jplev : level to treat if your file is 3D (spatial), has no influence if 
!!         your file is 2D !
!!         jplev > 0 = level to treat (ex : jplev = 1 will interpolate only 
!!                     surface field corresponding to the 1st level )
!!       ------------------------------------------------------------------
!!       | jplev = 0 : perform 3D interpolation (if input file is 3D) !!! |
!!       ------------------------------------------------------------------
!!
!! cf_x_in   : file containing the input grid (usually = cf_in) [char] 
!! cv_lon_in : name of longitude in the input grid file [char] 
!! cv_lat_in : name of latitude in the input grid file [char] 
!! 
!! cf_lsm_in : file containing the input land-sea mask [char]
!!             or specify 'missing_value' if a 'missing_value' netcdf
!!             attribute defines the mask on the input data field
!!             (not needed if "ldrown = .FALSE." --> '')
!!
!! cv_lsm_in : name of land-sea mask variable [char]
!!             (not needed if "ldrown = .FALSE." 
!!              or if cf_lsm_in = 'missing_value'--> '')
!!             by default ocean is flagged with value 1
!!             and continents are flagged with value 0
!!
!! ldrown : whether we call DROWN land filling procedure [logical] 
!!
!! ewper : east-west periodicity on the input file/grid [integer] 
!!         = -1 --> no periodicity 
!!         >= 0 --> periodicity with overlap of ewper points 
!! 
!! vmax : upper bound not to exceed for treated variable [real] 
!! vmin : lower bound not to exceed for treated variable [real] 
!!-------------------------------------------------------------------------- 
!! 
&ninput
ivect     = 0
lregin    = T
cf_in     = '/data/toto/my_file.nc'
cv_in     = 'temp'
cv_t_in   = 'time' 
jt1       = 0
jt2       = 0
jplev     = 1
cf_x_in   = '/data/toto/my_file.nc'
cv_lon_in = 'lon'
cv_lat_in = 'lat'
cf_lsm_in = 'missing_value'
cv_lsm_in = ''
ldrown    = T
ewper     = 0
vmax      =  1.E6
vmin      = -1.E6
/
!!
!!
!!
!!
!! ***********************************
!!  IF 3D INTERPOLATION ( jplev = 0 )
!! ***********************************
!!
!! Only mind if you do want to perform a 3D (spatial) interpolation
!! 
!! Mind only if you do want to perform a 3D interpolation !
!! First, make sure that jplev is set to 0 !
!!
!! cf_z_in  : file containing the input depth vector (associates a depth to a 
!!            given level). In most cases should be the same file than cf_x_in.
!! cv_z_in  : name of the variable for the input depth vector
!!
!! cf_z_out : file containing the output depth vector (associates a depth to a 
!!            given level). In most cases should be the same file than cf_x_in.
!! cv_z_out : name of the variable for the output depth vector in file 'cf_z_out'
!! cv_z_out_name: name you wish to give to depth variable in file to be created...
!!
!!
&n3d
cf_z_in  = '/data/toto/my_file.nc'
cv_z_in  = 'level'
cf_z_out = '/data/toto/mesh_zgr.nc'
cv_z_out = 'nav_lev'
cv_z_out_name = 'depth'
/
!! 
!!
!!
!!
!!
!! ***************************** 
!! Output Grid characteristics : 
!! ***************************** 
!! 
!! lregout : is the target grid regular ? [logical] 
!!           (ie : are output longitude and latitude 1D?) 
!!
!! cf_x_out   : file containing the target grid [char] 
!! cv_lon_out : name of longitude variable [char] 
!! cv_lat_out : name of latitude variable [char] 
!!
!! TRICK:  for interpolating onto a global regular spherical grid
!! ------  with a resolution of dx deg. of longitude and dy deg. of latitude
!!         * cf_x_out   = 'spheric' ! tells SOSIE to build a spherical output grid
!!         * cv_lon_out = '1.0'  ! your dx, here 1.0 deg. 
!!         * cv_lat_out = '1.0'  ! your dy, here 1.0 deg.
!!
!!
!! cf_lsm_out : file containing output land-sea mask [char] 
!!              MUST BE 3D for 3D interpolation!
!!              or specify 'missing_value' if a 'missing_value' netcdf
!!              attribute defines the mask on a field 'X' in file 'cf_x_out'
!!              (not needed if "lmout = .FALSE." --> '') 
!!
!! cv_lsm_out : name of land-sea mask variable in 'cf_lsm_out'     [char]
!!              or name of field 'X' in 'cf_x_out' if you specified 
!!              cf_lsm_out = 'missing_value'
!!              (not needed if "lmout = .FALSE." --> '') 
!!
!! lmout : whether to mask the interpolated field on the output file [logical] 
!!
!! rmaskvalue : missing value given to output field (for continents) [logical] 
!!
!! lct   : whether to control or not time variable [logical] 
!!         TRUE -> specify time array with starting time 't0' and step 't_stp' 
!!         FALSE -> same time array as in input file is used 
!! t0    : time to start (if lct is set to .TRUE.) [real] 
!! t_stp : time step (if lct is set to .TRUE.) [real] 
!!
!!
&noutput
lregout    = F
cf_x_out   = '/data/toto/coordinates_ORCA_R2.nc'
cv_lon_out = 'nav_lon'
cv_lat_out = 'nav_lat'
cf_lsm_out = '/data/toto/mask_ORCA2.nc'
cv_lsm_out = 'tmask'
lmout      = T
rmaskvalue = -9999.
lct        = F
t0         = 0.
t_stp      = 0.25
/
!! 
!! 
!! 
!! 
!! ******************************* 
!! Netcdf output characteristics : 
!! ******************************* 
!! 
!! This mostly deals with how the output file to be created is going to look like!
!!
!! cmethod  : the 2D interpolation method to be used
!!            use 'akima' if your input domain is regular (non-distorted grid)
!!            use 'bilin' otherwise, which is bilinear 2D interpolation
!!
!! *** Into the netcdf file to be created : *** 
!! cv_l_out : name for longitude on the output file [char] 
!! cv_p_out : name for latitude on the output file [char] 
!! cv_t_out : name for time on the output file [char] 
!! cv_out   : name for treated variable on the output file [char] 
!! cu_out   : treated variable units [char] 
!! cu_t     : time unit [char]
!! cln_out  : treated variable long name [char] 
!! cd_out   : directory to create output file to [char]
!! 
!! *** Name of the output file : *** 
!! csource  : short string to describe origin of the file [char] 
!! ctarget  : short string to describe the target grid [char] 
!! cextra   : short extra indication about the file [char] 
!! 
&nnetcdf
cmethod  = 'akima'
cv_l_out = 'nav_lon'
cv_p_out = 'nav_lat'
cv_t_out = 'time_counter'
cv_out   = 'temp'
cu_out   = 'deg. C'
cu_t     = 'unknown'
cln_out  = 'Temperature'
cd_out   = '.'
csource  = 'ERA40'
ctarget  = 'ORCA2'
cextra   = 'monthly_1998'
/ 
!!
         

A few examples

Into the "examples" sub-directory you will find namelist and data to test SOSIE on working configurations.
For each example you will find a commented and working namelist (namelist.exampleX) explaining the relevant namelist tuning.
We encourage you to use a software like ncview to visualize files to be interpolated and interpolated files.
The common approach to test a given example (# X):

        >> cd examples/
        >> tar zxvf data.tar.gz
        >> sosie.x -f namelist.exampleX
        
Example #1: basic 2D scalar field interpolation
Interpolation of Reynolds (2002) Long Term Mean SST onto the ORCA1 grid (illustrated on Fig. 3-4).
Do the interpolation:

      >> sosie.x -f namelist.example1
        
Check the newly created SST_360x180-ORCA1_REYNOLDS_LTM.nc


Example #2: 3D scalar field interpolation
3D interpolation of Levitus (1998) temperature climatology onto the ORCA1 grid (only march).

      >> sosie.x -f namelist.example2
        
Check the newly created temp_360x180-ORCA1_march.nc


Example #3: Interpolation from an irregular (ORCA1) to a regular lat-lon 1x1 deg. grid
2D interpolation of a SST snapshot from a random NEMO-ORCA1 simulation onto lat-lon 1x1 deg. grid using the bilinear algorithm (illustrated on Fig. 5-6).

      >> sosie.x -f namelist.example3
        
Check the newly created sst_ORCA1-1x1_test.nc


Example #4: Interpolation and correction of a 2D vector field from a regular lat-lon 1x1 deg. grid to an irregular grid (ORCA1)
As the ORCA family of grids gets distorted in the northern hemisphere it is necessary to correct (i.e. rotate) both components of the vector. In this example the input vector field is the wind at 10 from a few 6-hourly snapshots of the ERA-INTERIM re-analysis.

Do the "raw" interpolation for the zonal component of the wind:

      >> sosie.x -f namelist.example4_x
        
Do the "raw" interpolation for the meridional component of the wind:

      >> sosie.x -f namelist.example4_y
        
Now that uraw_1x1-deg-ORCA1_vct.nc and vraw_1x1-deg-ORCA1_vct.nc are created, time to correct:

 >> corr_vect.x -x u10 -y v10 -f namelist.example4_x -m data/mesh_mask_ORCA1_light.nc
        
Check the newly created u10_1x1-deg-ORCA1_vct.nc and u10_1x1-deg-ORCA1_vct.nc Example #5: Interpolation of a vector field from an irregular (ORCA) grid to a regular grid
Say that you want to interpolate the current field from an NEMO output to a regular grid.
  • file ORCA1_my_run_grid_U.nc contains the zonal current vozocrtx
  • file ORCA1_my_run_grid_V.nc contains the meridional current vomecrty
  • 
     >> corr_vect.x -I -x vozocrtx -y vomecrty -i ORCA1_my_run_grid_U.nc ORCA1_my_run_grid_V.nc \
                    -m mesh_mask_ORCA1_light.nc -t time_counter
            
    files ORCA1_my_run_grid_U_unrotated.nc and ORCA1_my_run_grid_V_unrotated.nc are generated and contain the "unrotated" version of each component of vector field, these fields can be normally interpolated to a regular grid using the bilinear method...





    Last edited by Laurent Brodeau, 2014/12/12