diff --git a/NEWS.txt b/NEWS.txt
index a7de360296d14befaf9d64dba821369ca0302ef4..190511a2c6f4015faaacac863f4d2857254938e7 100644
--- a/NEWS.txt
+++ b/NEWS.txt
@@ -1,6 +1,9 @@
 Latest changes
 ==============
 
+TELEMAC-2D/TELEMAC-3D/TOMAWAC/GAIA: add Nearest Neighbor search method for Time
+Series files
+
 Python: fix a bug with extract_contour on Windows occurring with some meshes
 and resulting in incorrect contour extraction.
 
diff --git a/documentation/telemac2d/user/latex/chap3.tex b/documentation/telemac2d/user/latex/chap3.tex
index 105740ceeb894a2d073817b24a344e8bbf7de443..e0236d6421bc9dc9ffe1b044bd84ba07c28fc9b5 100644
--- a/documentation/telemac2d/user/latex/chap3.tex
+++ b/documentation/telemac2d/user/latex/chap3.tex
@@ -527,8 +527,8 @@ written for certain points in the domain of the simulation.
 This feature has been introduced since release 8.5.
 In this way, one can obtain results with a high time frequency, without using
 a large amount of disk space.
-The time series file module interpolates data (using linear interpolation
-even if using quasi-bubble or quadratic elements)
+The time series file module interpolates data (using linear interpolation,
+by default, even if using quasi-bubble or quadratic elements)
 to the required output coordinates.
 \begin{itemize}
 \item \telkey{TIME SERIES FILE 1} The name of the output file to write.
@@ -599,6 +599,15 @@ Description of the file
 4th-6th line: output points; x and y coordinates, unique ID, and station name
 \end{verbatim}
 
+The keyword \telkey{INTERPOLATION METHOD FOR TIME SERIES} will define the type
+of interpolation to be used for extracting time series data, at the locations
+defined in the time series coordinates file (1 or 2).
+The two possibilities are:
+\begin{itemize}
+\item 0: the nearest neighbor search (that will also work for points outside the
+domain),
+\item 1: linear interpolation inside the triangle (default value).
+\end{itemize}
 
 \subsection{The listing printout}
 
diff --git a/documentation/telemac3d/user/latex/input_output.tex b/documentation/telemac3d/user/latex/input_output.tex
index 0f02b4f0772bcf8dcbc72aaabf50870a3e4cdfc0..82b8cd8052571af8014a3ee4b2cae98b2b21f253 100644
--- a/documentation/telemac3d/user/latex/input_output.tex
+++ b/documentation/telemac3d/user/latex/input_output.tex
@@ -750,8 +750,8 @@ written for certain points in the domain of the simulation.
 This feature has been introduced since release 8.5.
 In this way, one can obtain results with a high time frequency, without using
 a large amount of disk space.
-The time series file module interpolates data (using linear interpolation)
-to the required output coordinates.
+The time series file module interpolates data (using by default linear
+interpolation) to the required output coordinates.
 \begin{itemize}
 \item \telkey{2D TIME SERIES FILE} The name of the output file to write
 2D results.
@@ -821,6 +821,16 @@ Description of the file
 4th-6th line: output points; x and y coordinates, unique ID, and station name
 \end{verbatim}
 
+The keyword \telkey{INTERPOLATION METHOD FOR TIME SERIES} will define the type
+of interpolation to be used for extracting time series data, at the locations
+defined in the time series coordinates file.
+The two possibilities are:
+\begin{itemize}
+\item 0: the nearest neighbor search (that will also work for points outside the
+domain),
+\item 1: linear interpolation inside the triangle (default value).
+\end{itemize}
+
 \section{Results in longitude-latitude}
 If the geometry file is given in longitude-latitude, the result files
 can also be written with the coordinates in longitude-latitude
diff --git a/documentation/tomawac/user/latex/io.tex b/documentation/tomawac/user/latex/io.tex
index 21b494d82e9ead0ed8ed7b69098a8f47ae0df20f..0fc883cfeed8f8c54a677e8f15e315f142f1de1b 100644
--- a/documentation/tomawac/user/latex/io.tex
+++ b/documentation/tomawac/user/latex/io.tex
@@ -458,8 +458,8 @@ are written for certain points in the domain of the simulation.
 This feature has been introduced since release 8.5.
 In this way, one can obtain results with a high time frequency, without using
 a large amount of disk space.
-The time series file module interpolates data (using linear interpolation
-even if using quasi-bubble or quadratic elements)
+The time series file module interpolates data (using linear interpolation,
+by default, even if using quasi-bubble or quadratic elements)
 to the required output coordinates.
 \begin{itemize}
 \item \telkey{TIME SERIES FILE} The name of the output file to write.
@@ -523,6 +523,15 @@ Description of the file
 4th-6th line: output points; x and y coordinates, unique ID, and station name
 \end{verbatim}
 
+The keyword \telkey{INTERPOLATION METHOD FOR TIME SERIES} will define the type
+of interpolation to be used for extracting time series data, at the locations
+defined in the time series coordinates file.
+The two possibilities are:
+\begin{itemize}
+\item 0: the nearest neighbor search (that will also work for points outside the
+domain),
+\item 1: linear interpolation inside the triangle (default value).
+\end{itemize}
 
 \subsection{ The punctual or spectrum results file}
 \label{se:SpeFile}
diff --git a/examples/tomawac/3Dcoupling/t3d_3Dcoupling_NN.cas b/examples/tomawac/3Dcoupling/t3d_3Dcoupling_NN.cas
new file mode 100644
index 0000000000000000000000000000000000000000..3efcc6a6999e8dd7e7f8ae60945a72d3a9a740a3
--- /dev/null
+++ b/examples/tomawac/3Dcoupling/t3d_3Dcoupling_NN.cas
@@ -0,0 +1,67 @@
+TITLE = 'TELEMAC 3D : COUPLAGE TOMAWAC'
+/
+/----------------------------------------------------------------------/
+/ COUPLING WITH TOMAWAC                                                /
+/----------------------------------------------------------------------/
+FORTRAN FILE : user_fortran
+WAVE DRIVEN CURRENTS = YES
+COUPLING WITH = 'TOMAWACT3D'
+COUPLING PERIOD FOR TOMAWAC = 10
+TOMAWAC STEERING FILE = 'tom_3Dcouplittoral_NN.cas'
+BOUNDARY CONDITIONS FILE : geo_t3d_littoral.cli
+GEOMETRY FILE          : geo_littoral.slf
+3D RESULT FILE : r3d_littoralcoup.slf
+2D RESULT FILE : r2d_littoralcoup.slf
+
+TIME SERIES COORDINATES FILE : coordinate_t3d.txt
+2D TIME SERIES FILE : r2d_hist.slf
+3D TIME SERIES FILE : r3d_hist.slf
+INTERPOLATION METHOD FOR TIME SERIES = 0
+/
+/----------------------------------------------------------------------/
+/     OPTIONS GENERALES
+/----------------------------------------------------------------------/
+NUMBER OF TIME STEPS : 200
+TIME STEP            : 5.
+LISTING PRINTOUT PERIOD : 40
+GRAPHIC PRINTOUT PERIOD : 10
+VARIABLES FOR 2D GRAPHIC PRINTOUTS : U,V,S,H
+VARIABLES FOR 3D GRAPHIC PRINTOUTS : Z,U,V,W,US,VS,WS
+NUMBER OF HORIZONTAL LEVELS = 5
+INITIAL CONDITIONS : 'CONSTANT ELEVATION'
+INITIAL ELEVATION : 0.
+PRESCRIBED ELEVATIONS : 0.
+VELOCITY VERTICAL PROFILES = 2;2
+MASS-BALANCE = YES
+TIDAL FLATS = YES
+NON-HYDROSTATIC VERSION = NO
+/-------------------------
+/  CONVECTION-DIFFUSION-PROPAGATION
+/-------------------------
+/
+SCHEME FOR ADVECTION OF VELOCITIES = 1
+MAXIMUM NUMBER OF ITERATIONS FOR DIFFUSION OF VELOCITIES = 200
+ACCURACY FOR DIFFUSION OF VELOCITIES = 1.D-9
+/
+SOLVER FOR PROPAGATION = 1
+ACCURACY FOR PROPAGATION = 1.D-9
+/ diminuer si pas stable, jusqu'à 0
+FREE SURFACE GRADIENT COMPATIBILITY = 0.9
+/
+MASS-LUMPING FOR DEPTH = 1.
+IMPLICITATION FOR DEPTH = 1.
+IMPLICITATION FOR VELOCITIES = 1.
+MASS-LUMPING FOR DIFFUSION = 1.
+IMPLICITATION FOR DIFFUSION = 1.
+/
+HORIZONTAL TURBULENCE MODEL = 4
+VERTICAL TURBULENCE MODEL = 1
+
+/----------------------------/
+/     SOURCE TERMS
+/----------------------------/
+LAW OF BOTTOM FRICTION = 5
+FRICTION COEFFICIENT FOR THE BOTTOM = 0.001D0
+/
+/ DEFAULT VALUE UNTIL V8P5 KEPT FOR NON REGRESSION
+TREATMENT OF NEGATIVE DEPTHS = 1
diff --git a/examples/tomawac/3Dcoupling/tom_3Dcouplittoral_NN.cas b/examples/tomawac/3Dcoupling/tom_3Dcouplittoral_NN.cas
new file mode 100644
index 0000000000000000000000000000000000000000..74465cd42b315712ccf2827da5aae4b6353c40b5
--- /dev/null
+++ b/examples/tomawac/3Dcoupling/tom_3Dcouplittoral_NN.cas
@@ -0,0 +1,66 @@
+/DEBUGGER = 1
+/********************************************************************
+/             STEERING FILE
+/             TOMAWAC
+/
+/********************************************************************
+/--------------------------------------------------------------------
+/     FILES
+/--------------------------------------------------------------------
+LIMIT SPECTRUM MODIFIED BY USER = YES 
+GEOMETRY FILE = 'geo_littoral.slf'
+BOUNDARY CONDITIONS FILE = 'geo_tom_littoral.cli'
+2D RESULTS FILE = 'tom_couplittoral.slf'
+TIME SERIES COORDINATES FILE : coordinate_t3d.txt
+TIME SERIES FILE = tom_hist.slf
+INTERPOLATION METHOD FOR TIME SERIES = 0
+/--------------------------------------------------------------------
+/ INPUTS - OUTPUTS
+/--------------------------------------------------------------------
+TITLE = 'LITTORAL'
+PERIOD FOR GRAPHIC PRINTOUTS = 10
+VARIABLES FOR 2D GRAPHIC PRINTOUTS =
+HM0,DMOY,TPR5,ZF,WD,FX,FY,UX,UY
+PERIOD FOR LISTING PRINTOUTS = 20
+/--------------------------------------------------------------------
+/ DISCRETISATION
+/--------------------------------------------------------------------
+MINIMAL FREQUENCY = 0.05
+FREQUENTIAL RATIO = 1.1007
+NUMBER OF FREQUENCIES = 12
+NUMBER OF DIRECTIONS = 12
+TIME STEP = 10.
+NUMBER OF TIME STEP = 100
+/--------------------------------------------------------------------
+/ INITIAL CONDITIONS
+/--------------------------------------------------------------------
+INITIAL STILL WATER LEVEL = 0.
+TYPE OF INITIAL DIRECTIONAL SPECTRUM = 6
+INITIAL SIGNIFICANT WAVE HEIGHT = 1.
+INITIAL PEAK FREQUENCY = 0.125
+INITIAL PEAK FACTOR = 3.3
+INITIAL ANGULAR DISTRIBUTION FUNCTION = 1
+INITIAL WEIGHTING FACTOR FOR ADF = 1.0
+INITIAL MAIN DIRECTION 1 = 30.
+INITIAL DIRECTIONAL SPREAD 1 = 3
+/--------------------------------------------------------------------
+/ BOUNDARY CONDITIONS
+/--------------------------------------------------------------------
+TYPE OF BOUNDARY DIRECTIONAL SPECTRUM = 6
+BOUNDARY SIGNIFICANT WAVE HEIGHT = 1.
+BOUNDARY PEAK FREQUENCY = 0.125
+BOUNDARY PEAK FACTOR = 3.3
+BOUNDARY MAIN DIRECTION 1 = 30.
+BOUNDARY DIRECTIONAL SPREAD 1 = 3.
+/--------------------------------------------------------------------
+/ OPTIONS
+/--------------------------------------------------------------------
+MINIMUM WATER DEPTH = 0.05
+INFINITE DEPTH = NO
+CONSIDERATION OF SOURCE TERMS = YES
+CONSIDERATION OF A WIND = NO
+BOTTOM FRICTION DISSIPATION = 1
+/===========DEFERLEMENT
+NUMBER OF BREAKING TIME STEPS = 5
+/=========== 1 : Battjes et Janssen (1978)
+DEPTH-INDUCED BREAKING DISSIPATION = 1
diff --git a/examples/tomawac/3Dcoupling/vnv_3d_coupling.py b/examples/tomawac/3Dcoupling/vnv_3d_coupling.py
index 8552dddf0ae4dcf8c6f8ce4edac4fd84ccc28568..167c23f5b8214a26ee355c321d45f047d3c449ab 100644
--- a/examples/tomawac/3Dcoupling/vnv_3d_coupling.py
+++ b/examples/tomawac/3Dcoupling/vnv_3d_coupling.py
@@ -119,6 +119,20 @@ class VnvStudy(AbstractVnvStudy):
                        cas=cas)
         del cas
 
+        # littoral3D T3D+TOM scalar TEL2TOM + 3D coupling
+        self.add_study('vnv_13',
+                       'telemac3d',
+                       't3d_3Dcoupling_NN.cas')
+
+        # littoral3D T3D+TOM parallel TEL2TOM  + 3D coupling
+        cas = TelemacCas('t3d_3Dcoupling_NN.cas', get_dico('telemac3d'))
+        cas.set('PARALLEL PROCESSORS', 4)
+
+        self.add_study('vnv_14',
+                       'telemac3d',
+                       't3d_3Dcoupling_NN_par.cas',
+                       cas=cas)
+        del cas
 
     def _check_results(self):
         """
@@ -198,6 +212,22 @@ class VnvStudy(AbstractVnvStudy):
                             'vnv_4:T3DHI3',
                             eps=[1e-8])
 
+        self.check_epsilons('vnv_13:T3DHI2',
+                            'f2d_histcoup.slf',
+                            eps=[1e-10])
+
+        self.check_epsilons('vnv_13:T3DHI3',
+                            'f3d_histcoup.slf',
+                            eps=[1e-9])
+
+        self.check_epsilons('vnv_14:T3DHI2',
+                            'f2d_histcoup.slf',
+                            eps=[1e-10])
+
+        self.check_epsilons('vnv_14:T3DHI3',
+                            'f3d_histcoup.slf',
+                            eps=[1e-9])
+
         # TEL2TOM COUPLING Comparison with a reference for T3D.
         self.check_epsilons('vnv_5:T3DRES',
                             'f3d_littoral_diff.slf',
@@ -278,6 +308,16 @@ class VnvStudy(AbstractVnvStudy):
                             'vnv_12:WACRES',
                             eps=[1e-8])
 
+        # TEL2TOM COUPLING Comparison between sequential and parallel for T3D.
+        self.check_epsilons('vnv_13:T3DHI3',
+                            'vnv_14:T3DHI3',
+                            eps=[1e-15])
+
+        # TEL2TOM COUPLING Comparison between sequential and parallel for WAC.
+        self.check_epsilons('vnv_13:WACHI2',
+                            'vnv_14:WACHI2',
+                            eps=[1e-15])
+
     def _post(self):
         """
         Post-treatment processes
diff --git a/sources/api/api.cmdf b/sources/api/api.cmdf
index 1935e85f15dece93edb90acc838f664ffd15698d..9e0dd89df4cb056cd3a81665d0976baf4b8fc498 100644
--- a/sources/api/api.cmdf
+++ b/sources/api/api.cmdf
@@ -545,6 +545,7 @@ files: bief_def.f
   bief_allvec_in_block.f
   config_code.f
   outbief.f
+  nearest_node.F90
   out_history.F90
   gregtim.f
   jultim.f
diff --git a/sources/gaia/declarations_gaia.f b/sources/gaia/declarations_gaia.f
index 2389817d23586889a7035ead1aa6241deb073c56..9e187902ab919d0eec65b664187c80b0db9450d4 100644
--- a/sources/gaia/declarations_gaia.f
+++ b/sources/gaia/declarations_gaia.f
@@ -953,6 +953,10 @@
 !
       INTEGER, TARGET :: RESTART_RECORD
 !
+!>    INTERPOLATION METHOD FOR TIME SERIES
+!
+      INTEGER INTERPM
+!
 !-----------------------------------------------------------------------
 !
 !       5) LOGICAL VALUES
diff --git a/sources/gaia/gaia.dico b/sources/gaia/gaia.dico
index 5306eaf52885bdb9463d3f99846aa0bba5be3779..437e04eb3ab75abcbf856724530ca379bc11b65f 100644
--- a/sources/gaia/gaia.dico
+++ b/sources/gaia/gaia.dico
@@ -24,7 +24,7 @@
 /---SYSTEME-TELEMAC--------------------------------------GAIA-----------
 /
 / INDICES IN USE SO FAR:
-/ INTEGER INDEX USED: 1-3,7-9,11-12,17-20,23-30,33-34,39-47,51-52,54-54,58-67,251 OUT OF  45
+/ INTEGER INDEX USED: 1-3,7-9,11-12,17-20,23-30,33-34,39-47,51-54,58-67,251 OUT OF  46
 / REAL INDEX USED: 3-13,16-16,21-22,24-29,32-38,40-45,51-54,258-262 OUT OF  42
 / LOGICAL INDEX USED: 1-3,5-11,15-16,18-18,22-22,25-27,29-32,34-34,59 OUT OF  23
 / STRING INDEX USED: 1-9,11-13,16-16,22-22,24-27,30-31,33-42,53-53,55-55,59-60,64-64,69-70,100 OUT OF  38
@@ -1081,6 +1081,33 @@ AIDE1 =
 \item MED     : MED double precision format based on HDF5.
 \end{itemize}'
 /
+NOM = 'METHODE D''INTERPOLATION DES SERIES TEMPORELLES'
+NOM1 = 'INTERPOLATION METHOD FOR TIME SERIES'
+TYPE = INTEGER
+INDEX = 53
+TAILLE = 1
+DEFAUT = 1
+DEFAUT1 = 1
+MNEMO = 'INTERPM'
+CHOIX =
+'0="Plus proche voisin"';
+'1="Interpolation lineaire"'
+CHOIX1 =
+'0="Nearest neighbor search"';
+'1="Linear interpolation"'
+RUBRIQUE =
+'RESULTATS';'';''
+RUBRIQUE1 =
+'RESULTS';'';''
+NIVEAU = 0
+AIDE =
+'Methode d interpolation utilisee pour extraire les series temporelles
+aux coordonnees renseignees dans le
+\telkey{FICHIER DE COORDONNEES DES SERIES TEMPORELLES}.'
+AIDE1 =
+'Interpolation method used for the time series extraction at the
+locations provided in the \telkey{TIME SERIES COORDINATES FILE}.'
+/
 ////////////////////////////////////////////////////////////////////////
 /// 7-DATA FILES
 ////////////////////////////////////////////////////////////////////////
diff --git a/sources/gaia/gaia_init.F b/sources/gaia/gaia_init.F
index 615d95a57559db8d8963dec5238ceb264fc88be2..9e89f6da3933bd7885b72d7fafe39b3613abfa8c 100644
--- a/sources/gaia/gaia_init.F
+++ b/sources/gaia/gaia_init.F
@@ -791,9 +791,9 @@
 !
       CALL OUTHIST_PREPARE('GAI   ',GAI_FILES(GAICOO),
      &                     GAI_FILES(GAIHI2),MESH,NELEM,NPOIN,
-     &                     MARDAT,MARTIM,MAXVAR,SORLEO,TEXTE,VARSOR,1,
-     &                     PROTYP=PROTYP_TEL,LATIT=LATIT_TEL,
-     &                     LONGIT=LONGIT_TEL,
+     &                     MARDAT,MARTIM,MAXVAR,SORLEO,TEXTE,VARSOR,
+     &                     1,INTERPM,PROTYP=PROTYP_TEL,
+     &                     LATIT=LATIT_TEL,LONGIT=LONGIT_TEL,
      &                     KEEP_LONLAT=KEEP_LONLAT_TEL)
 
 !     PREPARES RESULTS
diff --git a/sources/gaia/lecdon_gaia.f b/sources/gaia/lecdon_gaia.f
index dd4206b4e375f2259c1ae62c431446191167f467..dee30f23133b70f116ed2b63b5ab50551614413c 100644
--- a/sources/gaia/lecdon_gaia.f
+++ b/sources/gaia/lecdon_gaia.f
@@ -707,6 +707,8 @@
       GAI_FILES(GAIHI2)%NAME=MOTCAR( ADRESS(4,25) )
       GAI_FILES(GAIHI2)%FMT = MOTCAR( ADRESS(4,26) )(1:8)
       CALL MAJUS(GAI_FILES(GAIHI2)%FMT)
+!     INTERPOLATION METHOD FOR THE TIME SERIES
+      INTERPM = MOTINT(ADRESS(1,53))
 
 !     RESULT FILE FORMAT FOR PREVIOUS SEDIMENTOLOGICAL
 !     COMPUTATION...
diff --git a/sources/telemac2d/declarations_telemac2d.f b/sources/telemac2d/declarations_telemac2d.f
index 23de91417215bac2016672c81558f2cf1dc22f9c..6965e104a4621fd9d48551ac96749f7299f1f99e 100644
--- a/sources/telemac2d/declarations_telemac2d.f
+++ b/sources/telemac2d/declarations_telemac2d.f
@@ -1558,6 +1558,10 @@
 !
       INTEGER OPTRTPF
 !
+!     INTERPOLATION METHOD FOR TIME SERIES
+!
+      INTEGER INTERPM
+!
 !-----------------------------------------------------------------------
 !
 !       5) LOGICAL VALUES
diff --git a/sources/telemac2d/lecdon_telemac2d.f b/sources/telemac2d/lecdon_telemac2d.f
index e3abf2eaaa7234ad3571414626e071d0dfcd324d..5c56a35971897083b3a159d5ad150e81bf9df8af 100644
--- a/sources/telemac2d/lecdon_telemac2d.f
+++ b/sources/telemac2d/lecdon_telemac2d.f
@@ -2097,6 +2097,8 @@
       CALL MAJUS(T2D_FILES(T2DHI1)%FMT)
       T2D_FILES(T2DHI2)%FMT = MOTCAR( ADRESS(4,109) )(1:8)
       CALL MAJUS(T2D_FILES(T2DHI2)%FMT)
+!     INTERPOLATION METHOD FOR THE TIME SERIES
+      INTERPM = MOTINT(ADRESS(1,122))
 !
 !    ===========================================================
 !
diff --git a/sources/telemac2d/telemac2d.cmdf b/sources/telemac2d/telemac2d.cmdf
index 238dedb5339bab7b4a64fe1efffee3d447b292da..58d520040769d943245d694cbc853f5305812e61 100644
--- a/sources/telemac2d/telemac2d.cmdf
+++ b/sources/telemac2d/telemac2d.cmdf
@@ -79,6 +79,7 @@ files: bief_def.f
   parcom2.f
   parcom2_seg.f
   parcom.f
+  nearest_node.F90
   out_history.F90
   gregtim.f
   jultim.f
diff --git a/sources/telemac2d/telemac2d.dico b/sources/telemac2d/telemac2d.dico
index c6d83e0f7a1b9c85458c07326646d56b11e2692a..2a8a4961b55f58d01e1fce7d2d6e980a917f6147 100644
--- a/sources/telemac2d/telemac2d.dico
+++ b/sources/telemac2d/telemac2d.dico
@@ -24,7 +24,7 @@
 /---SYSTEME-TELEMAC-------------------------------------------------------------
 /
 / INDICES IN USE SO FAR:
-/ INTEGER INDEX USED: 1-121 OUT OF 121
+/ INTEGER INDEX USED: 1-122 OUT OF 122
 / REAL INDEX USED: 1-98 OUT OF  98
 / LOGICAL INDEX USED: 3-20,22-40,46-51,53-64,66-68 OUT OF  58
 / STRING INDEX USED: 1-3,6-11,13-43,48-56,58-61,63-98,100-109 OUT OF  99
@@ -2021,6 +2021,34 @@ phrase "NORMAL END OF PROGRAM".
 In addition, the options \telkey{MASS-BALANCE} and
 \telkey{VALIDATION} are inhibited. Not recommended for use.'
 /
+NOM = 'METHODE D''INTERPOLATION DES SERIES TEMPORELLES'
+NOM1 = 'INTERPOLATION METHOD FOR TIME SERIES'
+TYPE = INTEGER
+INDEX = 122
+TAILLE = 1
+DEFAUT = 1
+DEFAUT1 = 1
+MNEMO = 'INTERPM'
+CHOIX =
+'0="Plus proche voisin"';
+'1="Interpolation lineaire"'
+CHOIX1 =
+'0="Nearest neighbor search"';
+'1="Linear interpolation"'
+RUBRIQUE =
+'ENVIRONNEMENT DE CALCUL';'INITIALISATION';'FICHIERS D ENTREE'
+RUBRIQUE1 =
+'COMPUTATION ENVIRONMENT';'INITIALIZATION';'INPUT FILES'
+NIVEAU = 0
+AIDE =
+'Methode d interpolation utilisee pour extraire les series temporelles
+aux coordonnees renseignees dans le
+\telkey{FICHIER DE COORDONNEES DES SERIES TEMPORELLES} (1 ou 2).'
+AIDE1 =
+'Interpolation method used for the time series extraction at the
+locations provided in the \telkey{TIME SERIES COORDINATES FILE}
+(1 or 2).'
+/
 NOM = 'VARIABLES A IMPRIMER'
 NOM1 = 'VARIABLES TO BE PRINTED'
 TYPE = STRING
diff --git a/sources/telemac2d/telemac2d_init.F b/sources/telemac2d/telemac2d_init.F
index 03cad01b37ab04f303f3b628597dc18522945366..893010455e7975e7c9d1aff81979b69df99d4d99 100644
--- a/sources/telemac2d/telemac2d_init.F
+++ b/sources/telemac2d/telemac2d_init.F
@@ -1163,15 +1163,15 @@
 
         CALL OUTHIST_PREPARE('T2D001',T2D_FILES(T2DCO1),
      &                       T2D_FILES(T2DHI1),MESH,NELEM,NPOIN,
-     &                       MARDAT,MARTIM,MAXVAR,SORLEO,TEXTE,VARSOR,1,
-     &                       PROTYP=PROTYP,LATIT=LAMBD0,LONGIT=PHI0,
-     &                       KEEP_LONLAT=KEEP_LONLAT)
+     &                       MARDAT,MARTIM,MAXVAR,SORLEO,TEXTE,VARSOR,
+     &                       1,INTERPM,PROTYP=PROTYP,LATIT=LAMBD0,
+     &                       LONGIT=PHI0,KEEP_LONLAT=KEEP_LONLAT)
 
         CALL OUTHIST_PREPARE('T2D002', T2D_FILES(T2DCO2),
      &                       T2D_FILES(T2DHI2),MESH,NELEM,NPOIN,
-     &                       MARDAT,MARTIM,MAXVAR,SORLEO,TEXTE,VARSOR,1,
-     &                       PROTYP=PROTYP,LATIT=LAMBD0,LONGIT=PHI0,
-     &                       KEEP_LONLAT=KEEP_LONLAT)
+     &                       MARDAT,MARTIM,MAXVAR,SORLEO,TEXTE,VARSOR,
+     &                       1,INTERPM,PROTYP=PROTYP,LATIT=LAMBD0,
+     &                       LONGIT=PHI0,KEEP_LONLAT=KEEP_LONLAT)
 !
 !=======================================================================
 !
diff --git a/sources/telemac3d/declarations_telemac3d.f b/sources/telemac3d/declarations_telemac3d.f
index 200f432f5286abc75a23d0cc8b4ae250d09bf4fb..aca7547c13508ac0fdafffdf15272ddcb2a31afa 100644
--- a/sources/telemac3d/declarations_telemac3d.f
+++ b/sources/telemac3d/declarations_telemac3d.f
@@ -1510,6 +1510,10 @@
 !
       INTEGER, TARGET :: T3DHI3
 !
+!     INTERPOLATION METHOD FOR TIME SERIES
+!
+      INTEGER INTERPM
+!
 !     FILES FOR DELWAQ
 !
       INTEGER T3DDL1,T3DDL2,T3DDL3,T3DDL4,T3DDL5,T3DDL6,T3DDL7,T3DDL8
diff --git a/sources/telemac3d/lecdon_telemac3d.F b/sources/telemac3d/lecdon_telemac3d.F
index abd7afb359436447c900de54c992cd72dabcabb6..89a3c3a4516b45e1acb1a66f41b5503e5db70b6b 100644
--- a/sources/telemac3d/lecdon_telemac3d.F
+++ b/sources/telemac3d/lecdon_telemac3d.F
@@ -1960,6 +1960,7 @@
 !     TIME SERIES FILES
       T3D_FILES(T3DHI2)%NAME=MOTCAR( ADRESS(4,37) )
       T3D_FILES(T3DHI3)%NAME=MOTCAR( ADRESS(4,38) )
+      INTERPM = MOTINT(ADRESS(1,78))
 !     OIL SPILL STEERING FILE
       T3D_FILES(T3DMIG)%NAME=MOTCAR( ADRESS(4,57) )
 !     FROM V9P0, IF AN OIL SPILL STEERING FILE IS DEFINED,
diff --git a/sources/telemac3d/telemac3d.cmdf b/sources/telemac3d/telemac3d.cmdf
index a061ac76224e28c1d7610b7485f5532ca655cfa0..e29506690a947424de263f0dbfab30208766ae1e 100644
--- a/sources/telemac3d/telemac3d.cmdf
+++ b/sources/telemac3d/telemac3d.cmdf
@@ -74,6 +74,7 @@ files: bief_def.f
   parcom2.f
   parcom2_seg.f
   parcom.f
+  nearest_node.F90
   out_history.F90
   gregtim.f
   jultim.f
diff --git a/sources/telemac3d/telemac3d.dico b/sources/telemac3d/telemac3d.dico
index 2f67cd5f5c0a583411b24aabf71ce1a7ce72471b..701cd773abb3d5b75a98c77a005243b636e6cf2f 100644
--- a/sources/telemac3d/telemac3d.dico
+++ b/sources/telemac3d/telemac3d.dico
@@ -23,7 +23,7 @@
 / AIDE1 : Help in English (LaTeX syntax)
 /---SYSTEME-TELEMAC-------------------------------------------------------------
 /
-/ INTEGER INDEX USED: 1-43,45-77,79-82,84-87,89-92,95-101,103-127 OUT OF 120
+/ INTEGER INDEX USED: 1-43,45-82,84-87,89-92,95-101,103-127 OUT OF 121
 / REAL INDEX USED: 1-42,51-107 OUT OF  99
 / LOGICAL INDEX USED: 2-2,5-18,21-28,31-34,51-51,58-61,68-72,76-76,81-84,86-87,90-93 OUT OF  48
 / STRING INDEX USED: 1-4,6-11,13-14,16-23,28-53,55-60,62-79,83-95,98-98,100-103 OUT OF  88
@@ -1565,6 +1565,33 @@ AIDE1 =
 \item MED     : MED double precision format based on HDF5.
 \end{itemize}'
 /
+NOM = 'METHODE D''INTERPOLATION DES SERIES TEMPORELLES'
+NOM1 = 'INTERPOLATION METHOD FOR TIME SERIES'
+TYPE = INTEGER
+INDEX = 78
+TAILLE = 1
+DEFAUT = 1
+DEFAUT1 = 1
+MNEMO = 'INTERPM'
+CHOIX =
+'0="Plus proche voisin"';
+'1="Interpolation lineaire"'
+CHOIX1 =
+'0="Nearest neighbor search"';
+'1="Linear interpolation"'
+RUBRIQUE =
+'ENVIRONNEMENT DE CALCUL';'SORTIE';'RESULTATS'
+RUBRIQUE1 =
+'COMPUTATION ENVIRONMENT';'OUTPUT';'RESULTS'
+NIVEAU = 0
+AIDE =
+'Methode d interpolation utilisee pour extraire les series temporelles
+aux coordonnees renseignees dans le
+\telkey{FICHIER DE COORDONNEES DES SERIES TEMPORELLES}.'
+AIDE1 =
+'Interpolation method used for the time series extraction at the
+locations provided in the \telkey{TIME SERIES COORDINATES FILE}.'
+/
 ////////////////////////////////////////////////////////////////////////
 /////// 1.3.2-LISTING
 ////////////////////////////////////////////////////////////////////////
diff --git a/sources/telemac3d/telemac3d_init.F b/sources/telemac3d/telemac3d_init.F
index 6ab72b121b50ebfef06aa539c45156e049e5b386..0914120eab52f2ebf2d5deb51ad9142afe49ee9c 100644
--- a/sources/telemac3d/telemac3d_init.F
+++ b/sources/telemac3d/telemac3d_init.F
@@ -1029,15 +1029,15 @@
 !     2D TIME SERIES FILE
       CALL OUTHIST_PREPARE('T3D2DH',T3D_FILES(T3DCOO),
      &                     T3D_FILES(T3DHI2),MESH2D,NELEM2,NPOIN2,
-     &                     MARDAT,MARTIM,MAXVAR,SORG2D,TEXTE,VARSOR,1,
-     &                     PROTYP=PROTYP,LATIT=LATIT,LONGIT=LONGIT,
-     &                     KEEP_LONLAT=KEEP_LONLAT)
+     &                     MARDAT,MARTIM,MAXVAR,SORG2D,TEXTE,VARSOR,
+     &                     1,INTERPM,PROTYP=PROTYP,LATIT=LATIT,
+     &                     LONGIT=LONGIT,KEEP_LONLAT=KEEP_LONLAT)
 !
 !     3D TIME SERIES FILE
       CALL OUTHIST_PREPARE('T3D3DH',T3D_FILES(T3DCOO),
      &                     T3D_FILES(T3DHI3),MESH2D,NELEM2,NPOIN2,
      &                     MARDAT,MARTIM,MAXVA3,SORG3D,TEXT3,VARSO3,
-     &                     NPLAN,PROTYP=PROTYP,LATIT=LATIT,
+     &                     NPLAN,INTERPM,PROTYP=PROTYP,LATIT=LATIT,
      &                     LONGIT=LONGIT,KEEP_LONLAT=KEEP_LONLAT)
 !
 !     PREPARES THE 2D AND 3D OUTPUT FOR INITIAL CONDITIONS
diff --git a/sources/tomawac/declarations_tomawac.f b/sources/tomawac/declarations_tomawac.f
index a388673ce15d1f66bb145eb5acbda95d667c0241..b8d1f4ecb47fa0a4f802327351cd5f96dfc4e1c2 100644
--- a/sources/tomawac/declarations_tomawac.f
+++ b/sources/tomawac/declarations_tomawac.f
@@ -1322,6 +1322,10 @@
 !
       INTEGER, TARGET :: WACHI2
 !
+!     INTERPOLATION METHOD FOR TIME SERIES
+!
+      INTEGER INTERPM
+!
 !     COORDINATES OF POINTS FOR TIME SERIES FILE
 !
       INTEGER, TARGET :: WACCOO
diff --git a/sources/tomawac/lecdon_tomawac.f b/sources/tomawac/lecdon_tomawac.f
index 83737ae0a1efa39541ab95e41e9ee1c6e70fade3..1df543f9a1852501b1179c9396baf22678d4f571 100644
--- a/sources/tomawac/lecdon_tomawac.f
+++ b/sources/tomawac/lecdon_tomawac.f
@@ -568,6 +568,8 @@
  !    FORMAT OF THE TIME SERIES FILE
       WAC_FILES(WACHI2)%FMT = MOTCAR( ADRESS(4,30) )(1:8)
       CALL MAJUS(WAC_FILES(WACHI2)%FMT)
+!     INTERPOLATION METHOD FOR THE TIME SERIES
+      INTERPM = MOTINT(ADRESS(1,58))
       IF( WAC_FILES(WACRBI)%NAME.NE.' ') THEN
 ! ONE WANTS TO HAVE A GLOBAL RESULT
         GLOB=.TRUE.
diff --git a/sources/tomawac/tomawac.cmdf b/sources/tomawac/tomawac.cmdf
index f772fa90ea36a42b04d8a0bcb0cf2c79098c6650..f002b5bf3ef779b58a7ab701c69248ab97be921a 100644
--- a/sources/tomawac/tomawac.cmdf
+++ b/sources/tomawac/tomawac.cmdf
@@ -47,6 +47,7 @@ files: bief_def.f
   os.f
   algae_transp.f
   deall_bief.f
+  nearest_node.F90
   out_history.F90
   allblo.f
   bief_allvec_in_block.f
diff --git a/sources/tomawac/tomawac.dico b/sources/tomawac/tomawac.dico
index bf29e0d6685919d4a85ec1c94e6655354a39ea30..860e9d357c25b5439e9eef9df02bb5017f771cb7 100644
--- a/sources/tomawac/tomawac.dico
+++ b/sources/tomawac/tomawac.dico
@@ -1,5 +1,5 @@
 &DYN
-/ INTEGER INDEX USED: 1-57 OUT OF  57
+/ INTEGER INDEX USED: 1-58 OUT OF  58
 / REAL INDEX USED: 1-100 OUT OF 100
 / LOGICAL INDEX USED: 1-2,5-25,97 OUT OF  24
 / STRING INDEX USED: 1-15,17-17,23-32,34-36,39-47,50-51,100-101 OUT OF  42
@@ -6303,7 +6303,34 @@ AIDE1 =
 \item SERAFIND: classical double precision format in \tel,
 \item MED     : MED double precision format based on HDF5.
 \end{itemize}'
-
+/
+NOM = 'METHODE D''INTERPOLATION DES SERIES TEMPORELLES'
+NOM1 = 'INTERPOLATION METHOD FOR TIME SERIES'
+TYPE = INTEGER
+INDEX = 58
+TAILLE = 1
+DEFAUT = 1
+DEFAUT1 = 1
+MNEMO = 'INTERPM'
+CHOIX =
+'0="Plus proche voisin"';
+'1="Interpolation lineaire"'
+CHOIX1 =
+'0="Nearest neighbor search"';
+'1="Linear interpolation"'
+RUBRIQUE =
+'ENVIRONNEMENT DE CALCUL';'SORTIE';'RESULTATS'
+RUBRIQUE1 =
+'COMPUTATION ENVIRONMENT';'OUTPUT';'RESULTS'
+NIVEAU = 0
+AIDE =
+'Methode d''interpolation utilisee pour extraire les series temporelles
+aux coordonnees renseignees dans le
+\telkey{FICHIER DE COORDONNEES DES SERIES TEMPORELLES}.'
+AIDE1 =
+'Interpolation method used for the time series extraction at the
+locations provided in the \telkey{TIME SERIES COORDINATES FILE}.'
+/
 NOM = 'UNITE DE TEMPS DU FICHIER DES SPECTRES IMPOSES'
 NOM1 = 'TIME UNIT OF IMPOSED SPECTRA FILE'
 TYPE = REAL
diff --git a/sources/tomawac/wac_init.F b/sources/tomawac/wac_init.F
index c9f8d04a18ae17449dda55bc998ad26553aa3baf..496f6c5383b713721ae0f0011f777b369deb8e23 100644
--- a/sources/tomawac/wac_init.F
+++ b/sources/tomawac/wac_init.F
@@ -470,8 +470,8 @@
 
       CALL OUTHIST_PREPARE('WAC   ',WAC_FILES(WACCOO),
      &                     WAC_FILES(WACHI2),MESH,NELEM2,NPOIN2,
-     &                     DATE,TIME,MAXVAR,SORLEO,TEXTE,VARSOR,1,
-     &                     SPHE=SPHE)
+     &                     DATE,TIME,MAXVAR,SORLEO,TEXTE,VARSOR,
+     &                     1,INTERPM,SPHE=SPHE)
 
       CALL HISTPERIOD('WAC   ',AT,IMPRES_HIS)
       IF(IMPRES.OR.IMPRES_HIS.OR.LT.EQ.0) THEN
diff --git a/sources/utils/bief/nearest_node.F90 b/sources/utils/bief/nearest_node.F90
new file mode 100644
index 0000000000000000000000000000000000000000..a62d1cc1d18ea3452e8b9f1dda9d646a0eafe535
--- /dev/null
+++ b/sources/utils/bief/nearest_node.F90
@@ -0,0 +1,110 @@
+
+!                   ***********************
+                    SUBROUTINE NEAREST_NODE &
+!                   ***********************
+!
+     (IP,XP,YP,X,Y,NP,NPOIN)
+!
+!***********************************************************************
+! BIEF
+!***********************************************************************
+!
+!>@brief    Inpired on PROXIM routine, 
+!           this works for points even outside of the domain
+!
+
+!history T SAILLOUR (EU JRC)
+!+        10/07/2024
+!+        v9.1
+!+   Add nearest point interpolation method
+!+   (Works also for points outsite the domain)
+!
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!>@param[in,out] IP     index of the nearest node in the mesh
+!>@param[in]     XP     x-coordinates of the point
+!>@param[in]     YP     y-coordinates of the point
+!>@param[in]     X      x-coordinates of the mesh
+!>@param[in]     Y      y-coordinates of the mesh
+!>@param[in]     NP     number of points to interpolate
+!>@param[in]     NPOIN  number of points in the mesh
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!
+      USE BIEF
+      USE DECLARATIONS_SPECIAL
+      USE INTERFACE_PARALLEL
+!
+      IMPLICIT NONE
+!
+!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
+!
+      INTEGER, INTENT(IN)    :: NP,NPOIN
+      INTEGER, INTENT(INOUT) :: IP(NP)
+!
+      DOUBLE PRECISION, INTENT(IN) :: XP(NP),YP(NP),X(NPOIN),Y(NPOIN)
+!
+!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
+!
+      INTEGER I,ICLOSE,K
+      DOUBLE PRECISION X1,Y1,DIST2,D2,XX,YY
+!
+      INTRINSIC SQRT
+!
+!-----------------------------------------------------------------------
+!
+      IF(NP.GT.0) THEN
+!
+!       LOOP OVER THE OBSERVATION POINTS
+!
+        DO K=1,NP
+          IP(K)=0
+          DIST2=HUGE(1.D0)
+          XX=-HUGE(1.D0)
+          YY=-HUGE(1.D0)
+!
+!         LOOP ON THE POINTS TO FIND THE NEAREST POINT IN THE (SUB-)DOMAIN:
+!
+          DO I=1,NPOIN
+!           TAKES THE NEAREST NODE
+            D2=(XP(K)-X(I))**2+(YP(K)-Y(I))**2
+            IF(D2.LT.DIST2) THEN
+              IP(K)=I
+              XX=X(IP(K))
+              YY=Y(IP(K))
+              DIST2=D2
+            ENDIF
+          ENDDO
+!
+!         IN PARALLEL SEVERAL NEAREST POINTS MAY HAVE BEEN FOUND, 
+!         IF A POINT IN THE GLOBAL DOMAIN IS CLOSER THAN IP(K)
+!         WE ASSIGN IP(K) = 0
+          X1=XX
+          Y1=YY
+          IF(NCSIZE.GT.1) THEN
+            IF(DIST2.EQ.P_MIN(DIST2)) THEN
+              X1=XX
+              Y1=YY
+            ELSE 
+              X1=-HUGE(1.D0)
+              Y1=-HUGE(1.D0)
+              IP(K) = 0
+            ENDIF
+            X1=P_MAX(X1)
+            Y1=P_MAX(Y1)
+          ELSE
+            X1=YY
+            Y1=YY
+          ENDIF
+          WRITE(LU,*)
+          WRITE(LU,*) 'SOURCE POINT ',K,'PUT ON POINT'
+          WRITE(LU,*) X1,' AND ',Y1
+          D2 = SQRT(P_MIN(DIST2))
+          WRITE(LU,*) 'LOCATED AT ',D2,' DEGREES or METERS'
+!         LINE FEED FOR THE LISTING
+          IF(K.EQ.NP) WRITE(LU,*)
+        ENDDO
+      ENDIF
+!
+!-----------------------------------------------------------------------
+!
+      RETURN
+      END SUBROUTINE NEAREST_NODE
diff --git a/sources/utils/bief/out_history.F90 b/sources/utils/bief/out_history.F90
index b3bd98b8f1ece029bd60c852607be0a822a7c39f..4bb23ed93d389abab5854cc23540c6bc233a00da 100644
--- a/sources/utils/bief/out_history.F90
+++ b/sources/utils/bief/out_history.F90
@@ -96,9 +96,9 @@
 !+
 !
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!>@param [in,out] LEO     Write history file at this time step
-!>@param [in]     TIME    Time of current time step
 !>@param [in]     MODNAME Name of the calling module
+!>@param [in]     TIME    Time of current time step
+!>@param [in,out] LEO     Write history file at this time step
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !
 !
@@ -365,7 +365,7 @@
 !                   **************************
 !
       (MODNAME,COORFILE,OUTFILE,MESH2D,NELEM2,NPOIN2,DATE,TIME,MAXVAR, &
-       SORG2D,TEXTE,VARSOR,NPLAN,PROTYP,LATIT,LONGIT,KEEP_LONLAT,SPHE)
+       SORG2D,TEXTE,VARSOR,NPLAN,INTERPM,PROTYP,LATIT,LONGIT,KEEP_LONLAT,SPHE)
 !
 !***********************************************************************
 ! BIEF
@@ -398,6 +398,7 @@
 !>@param[in]  TEXTE    Variable names and units
 !>@param[in]  VARSOR   Output data
 !>@param[in]  NPLAN    Number of vertical planes
+!>@param[in]  INTERPM  Interpolation method
 !>@param[in]  PROTYP   Spatial projection type (3: longitude/latitude)
 !>@param[in]  LATIT    Latitude of origin point (Mercator conversion)
 !>@param[in]  LONGIT   Longitude of origin point (Mercator conversion)
@@ -418,6 +419,7 @@
       INTEGER , INTENT(IN) :: NELEM2,NPOIN2,NPLAN
       INTEGER, INTENT(IN) :: MAXVAR
       INTEGER,DIMENSION(3), INTENT(IN) :: DATE,TIME
+      INTEGER, INTENT(IN) :: INTERPM
       LOGICAL, INTENT(IN) :: SORG2D(MAXVAR)
       CHARACTER(LEN=32),INTENT(IN) :: TEXTE(MAXVAR)
       TYPE(BIEF_OBJ), INTENT(IN) :: VARSOR
@@ -643,7 +645,7 @@
 
       CALL INTERP_POINT_PREPARE(TCDAT(IFILE)%XHISTOUT, &
            TCDAT(IFILE)%YHISTOUT,TCDAT(IFILE)%NRPOINT,INTDAT(IFILE), &
-           INTERP_METHOD,MESH2D,MODNAME,PROJECTION,SPHERI)
+           INTERPM,MESH2D,MODNAME,PROJECTION,SPHERI)
 
       ! Determine which points to consider
       INTDAT(IFILE)%MASKARR = INTDAT(IFILE)%POINTLIST(:,1).NE.0
@@ -1071,7 +1073,7 @@
 !>@param[in]      NPINT  Number of points to interpolate
 !>@param[in]      XINT   Structure with coordinates to use
 !>@param[in,out]  INTVAR Structure with data for interpolation
-!>@param[in]      INTERP_METHOD 1 if linear interpolation, 0 if nearest
+!>@param[in]      INTERPM 1 if linear interpolation, 0 if nearest
 !                               neighbour interpolation 
 !>@param[in]      MESH2D Structure with the 2d mesh
 !>@param [in]     MODNAME Name of the calling module
@@ -1105,76 +1107,79 @@
 !
 !-----------------------------------------------------------------------
 !
-      ! Preparxation
+      ! Preparation
       NELEM2 = MESH2D%NELEM
       IKLE2=>MESH2D%IKLE%I
       INTVAR%POINTLIST  = 0
       INTVAR%WEIGHTS = 0.D0
 
 ! Loop over all points that need to interpolated
-
-      DO I=1,NPINT
-
-        XP = XINT(I)
-        YP = YINT(I)
-        ! Loop over all elements in the mesh
-        DO J=1,NELEM2
-          J2 = J + NELEM2
-          J3 = J + 2*NELEM2
-
-          ! Get coordinates of the triangle
-          XTRI(1) = MESH2D%X%R(IKLE2(J))
-          XTRI(2) = MESH2D%X%R(IKLE2(J2))
-          XTRI(3) = MESH2D%X%R(IKLE2(J3))
-          YTRI(1) = MESH2D%Y%R(IKLE2(J))
-          YTRI(2) = MESH2D%Y%R(IKLE2(J2))
-          YTRI(3) = MESH2D%Y%R(IKLE2(J3))
-
-          ! Check if a point is in a trinagle
-          CALL IN_TRIANGLE(XP,YP,XTRI,YTRI,INTRI,MODNAME,PROTYP,SPHERI)
-
-          IF(INTRI) THEN
-            ! Store points
-            INTVAR%POINTLIST(I,1) = IKLE2(J)
-            INTVAR%POINTLIST(I,2) = IKLE2(J2)
-            INTVAR%POINTLIST(I,3) = IKLE2(J3)
-            ! Determine weights for linear interpolation
-            ! Area  p-2-3
-            XTRI(1) = XP
-            YTRI(1) = YP
-            CALL TRI_AREA(XTRI,YTRI,W1)
-            !Area p-1-3: point
-            XTRI(2) = MESH2D%X%R(IKLE2(J))
-            YTRI(2) = MESH2D%Y%R(IKLE2(J))
-            CALL TRI_AREA(XTRI,YTRI,W2)
-            !Area p-1-2: point
-            XTRI(3) = MESH2D%X%R(IKLE2(J2))
-            YTRI(3) = MESH2D%Y%R(IKLE2(J2))
-            CALL TRI_AREA(XTRI,YTRI,W3)
-            W = W1 + W2 + W3
-
-            IF(INTERP_METHOD.EQ.1) THEN
-!             Linear interpolation
+      IF(INTERP_METHOD.EQ.1) THEN
+        DO I=1,NPINT
+          XP = XINT(I)
+          YP = YINT(I)
+
+          ! Loop over all elements in the mesh
+          DO J=1,NELEM2
+            J2 = J + NELEM2
+            J3 = J + 2*NELEM2
+
+            ! Get coordinates of the triangle
+            XTRI(1) = MESH2D%X%R(IKLE2(J))
+            XTRI(2) = MESH2D%X%R(IKLE2(J2))
+            XTRI(3) = MESH2D%X%R(IKLE2(J3))
+            YTRI(1) = MESH2D%Y%R(IKLE2(J))
+            YTRI(2) = MESH2D%Y%R(IKLE2(J2))
+            YTRI(3) = MESH2D%Y%R(IKLE2(J3))
+
+            ! Check if a point is in a trinagle
+            CALL IN_TRIANGLE(XP,YP,XTRI,YTRI,INTRI,MODNAME,PROTYP,SPHERI)
+
+            IF(INTRI) THEN
+              ! Store points
+              INTVAR%POINTLIST(I,1) = IKLE2(J)
+              INTVAR%POINTLIST(I,2) = IKLE2(J2)
+              INTVAR%POINTLIST(I,3) = IKLE2(J3)
+              ! Determine weights for linear interpolation
+              ! Area  p-2-3
+              XTRI(1) = XP
+              YTRI(1) = YP
+              CALL TRI_AREA(XTRI,YTRI,W1)
+              ! Area p-1-3: point
+              XTRI(2) = MESH2D%X%R(IKLE2(J))
+              YTRI(2) = MESH2D%Y%R(IKLE2(J))
+              CALL TRI_AREA(XTRI,YTRI,W2)
+              ! Area p-1-2: point
+              XTRI(3) = MESH2D%X%R(IKLE2(J2))
+              YTRI(3) = MESH2D%Y%R(IKLE2(J2))
+              CALL TRI_AREA(XTRI,YTRI,W3)
+              W = W1 + W2 + W3
+              ! Linear interpolation
               INTVAR%WEIGHTS (I,1) = W1/W
               INTVAR%WEIGHTS (I,2) = W2/W
               INTVAR%WEIGHTS (I,3) = W3/W
-            ELSE
-              ! Nearest neighbour
-              IF((W1.GE.W2).AND.(W1.GE.W3)) THEN
-                INTVAR%WEIGHTS (I,1) = 1.D0
-              ENDIF
-              IF((W2.GE.W1).AND.(W2.GE.W3)) THEN
-                INTVAR%WEIGHTS (I,2) = 1.D0
-              ENDIF
-              IF((W3.GE.W2).AND.(W3.GE.W1)) THEN
-                INTVAR%WEIGHTS (I,3) = 1.D0
-              ENDIF
+              ! Continue to next point of the list
+              EXIT
             ENDIF
-            ! Continue to next point of the list
-            EXIT
-          ENDIF
+          ENDDO
         ENDDO
-      ENDDO
+      ELSEIF(INTERP_METHOD.EQ.0) THEN
+        ! Store the index of the nearest node
+        CALL NEAREST_NODE(INTVAR%POINTLIST,XINT,YINT,MESH2D%X%R, &
+                          MESH2D%Y%R,NPINT,MESH2D%NPOIN)
+        ! Set the weight of the nearest node to 1 and others to 0
+        DO I=1,NPINT
+          ! we need to assign I2 and I3 too, otherwise output
+          ! interpolation will not work
+          INTVAR%POINTLIST(I,2) = INTVAR%POINTLIST(I,1)
+          INTVAR%POINTLIST(I,3) = INTVAR%POINTLIST(I,1)
+          INTVAR%WEIGHTS(I,1) = 1.D0
+        ENDDO
+      ELSE
+        WRITE(LU,*) 'UNKNOWN INTERPOLATION METHOD FOR TIME SERIES'
+        WRITE(LU,*) 'USE 0 (NEAREST NEIGHBOR) OR 1 (LINEAR)'
+        CALL PLANTE(1)
+      ENDIF
 !
 !-----------------------------------------------------------------------
 !
@@ -1200,11 +1205,11 @@
 !+
 !
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!>@param[in]      DATA_IN  Data to interpolate
 !>@param[in]      INTVAR   Structure with data for interpolation
+!>@param[in]      NRPOINTS Size of the input data
+!>@param[in]      DATA_IN  Data to interpolate
 !>@param[in,out]  DATA_OUT Interpolated data
 !>@param[in]      OUT_SIZE Size of the output data
-!>@param[in]      NRPOINTS Size of the input data
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !
 !
@@ -1403,5 +1408,4 @@
 !
       RETURN
       END SUBROUTINE TRI_AREA
-
-      END MODULE OUT_HISTORY
+      END MODULE OUT_HISTORY
\ No newline at end of file