diff --git a/NEWS.txt b/NEWS.txt
index faeab047ff9388a4c444af84a35daddea84e055a..8e553e71aeff5f3957d4c12e5489e611572781ae 100644
--- a/NEWS.txt
+++ b/NEWS.txt
@@ -1,6 +1,8 @@
 Latest changes
 ==============
 
+TOMAWAC: fix time series for wave directions
+
 MASCARET: Reference Documentation has been added
 
 TELEMAC-2D: better time step managment in the finite volume kernel.
diff --git a/examples/tomawac/3Dcoupling/coordinate_tom.txt b/examples/tomawac/3Dcoupling/coordinate_tom.txt
new file mode 100644
index 0000000000000000000000000000000000000000..cfa3d5d87df08c9031f569aab989307bd4674d63
--- /dev/null
+++ b/examples/tomawac/3Dcoupling/coordinate_tom.txt
@@ -0,0 +1,22 @@
+1 8
+0 1000 50
+500    25       1 tom_01
+500    50       2 tom_02
+500    75       3 tom_03
+500   100       4 tom_04
+500   125       5 tom_05
+500   150       6 tom_06
+500   175       7 tom_07
+500   200       8 tom_08
+
+Description of the file
+1st line: number of periods and number of points
+2nd line: period 1: start time, end time and interval (in seconds)
+3rd-10th line: output points; x coordinate, y coordinate, station number, and station name
+
+The x and y coordinate are used for linear interpolation.
+The station number is an index, which is written to IKLE. (So you can easily know which station is in the output file)
+The station name is not used by TELEMAC, but is to increase the readibility of this file.
+The rest of the file are comments
+
+This file is specially made to test the use of two different output periods
diff --git a/examples/tomawac/3Dcoupling/fom_histcoup.slf b/examples/tomawac/3Dcoupling/fom_histcoup.slf
new file mode 100644
index 0000000000000000000000000000000000000000..fc4173b0a2a96acb3539c4ce05e90bd139f18fd3
--- /dev/null
+++ b/examples/tomawac/3Dcoupling/fom_histcoup.slf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c3236c59fb17aed9c2cb7ea871aa6c07dbb7d348d79df70626579e7d2c2ba5ef
+size 1068
diff --git a/examples/tomawac/3Dcoupling/tom_3Dcouplittoral.cas b/examples/tomawac/3Dcoupling/tom_3Dcouplittoral.cas
index 93f9102adc7244883989216eef3a8bb521b5bf77..979e8adff819c80ae9cbb03b30e0e22d478f3b33 100644
--- a/examples/tomawac/3Dcoupling/tom_3Dcouplittoral.cas
+++ b/examples/tomawac/3Dcoupling/tom_3Dcouplittoral.cas
@@ -11,6 +11,8 @@ 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_tom.txt'
+TIME SERIES FILE = 'tom_couplittoral_hist.slf'
 /--------------------------------------------------------------------
 / INPUTS - OUTPUTS
 /--------------------------------------------------------------------
diff --git a/examples/tomawac/3Dcoupling/vnv_3d_coupling.py b/examples/tomawac/3Dcoupling/vnv_3d_coupling.py
index 167c23f5b8214a26ee355c321d45f047d3c449ab..94e9ba7dbe84445a3a5c3378e71f0861f6a0e95e 100644
--- a/examples/tomawac/3Dcoupling/vnv_3d_coupling.py
+++ b/examples/tomawac/3Dcoupling/vnv_3d_coupling.py
@@ -211,6 +211,18 @@ class VnvStudy(AbstractVnvStudy):
         self.check_epsilons('vnv_3:T3DHI3',
                             'vnv_4:T3DHI3',
                             eps=[1e-8])
+                            
+        self.check_epsilons('vnv_3:WACHI2',
+                            'fom_histcoup.slf',
+                            eps=[1e-8])
+
+        self.check_epsilons('vnv_4:WACHI2',
+                            'fom_histcoup.slf',
+                            eps=[1e-8])
+
+        self.check_epsilons('vnv_3:WACHI2',
+                            'vnv_4:WACHI2',
+                            eps=[1e-8])
 
         self.check_epsilons('vnv_13:T3DHI2',
                             'f2d_histcoup.slf',
diff --git a/sources/gaia/deall_gaia.f b/sources/gaia/deall_gaia.f
index de02118806dc581e820fc172d164adedec33ef3b..7a8c30fe8d8942ea34093d5e509e76aaa2db38fc 100644
--- a/sources/gaia/deall_gaia.f
+++ b/sources/gaia/deall_gaia.f
@@ -46,6 +46,7 @@
 
 
       CALL BIEF_DEALLOBJ(VARSOR)
+      DEALLOCATE(ISDIR_GAIA)
 !     Deallocating all the blocks first
       CALL BIEF_DEALLOBJ(S     )
       CALL BIEF_DEALLOBJ(E     )
diff --git a/sources/gaia/declarations_gaia.f b/sources/gaia/declarations_gaia.f
index 9e187902ab919d0eec65b664187c80b0db9450d4..35eb291b1c9c2fbb322e35f9fea165231e4687ed 100644
--- a/sources/gaia/declarations_gaia.f
+++ b/sources/gaia/declarations_gaia.f
@@ -537,6 +537,10 @@
 !
       TYPE(BIEF_OBJ), TARGET :: VARSOR
 !
+!     Whether an output variable is a direction
+!
+      LOGICAL, ALLOCATABLE, TARGET :: ISDIR_GAIA(:)
+!
 !>    Vertical sorting profile: fraction for each layer, class, point
 !
       DOUBLE PRECISION,DIMENSION(:,:,:),TARGET,ALLOCATABLE::PRO_F
diff --git a/sources/gaia/gaia_init.F b/sources/gaia/gaia_init.F
index 9e89f6da3933bd7885b72d7fafe39b3613abfa8c..9d5678c29d422e1d5025cc7b68e993c67f7684d5 100644
--- a/sources/gaia/gaia_init.F
+++ b/sources/gaia/gaia_init.F
@@ -800,7 +800,7 @@
       CALL PREDES(0,AT0,.TRUE.,CODE,LISTCOUNT)
       IF(YAHIST) THEN
         CALL OUTHIST('GAI   ',GAI_FILES(GAIHI2),AT0,NPOIN,1,MAXVAR,
-     &               SORLEO,VARSOR)
+     &               SORLEO,VARSOR,ISDIR_GAIA)
       ENDIF
 !
 !     PRINTS OUT THE RESULTS
diff --git a/sources/gaia/gaia_write_results.f b/sources/gaia/gaia_write_results.f
index 01fcaba5db355984c387594021d427348b584223..4f6582a5282a38ee1817cbd7b2c4fcbdcdb9a816 100644
--- a/sources/gaia/gaia_write_results.f
+++ b/sources/gaia/gaia_write_results.f
@@ -49,7 +49,7 @@
       
       IF(YAHIST) THEN
         CALL OUTHIST('GAI   ',GAI_FILES(GAIHI2),T_TEL,NPOIN,1,MAXVAR,
-     &               SORLEO,VARSOR)
+     &               SORLEO,VARSOR,ISDIR_GAIA)
       ENDIF
            
       IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE BIEF_DESIMP'
diff --git a/sources/gaia/point_gaia.f b/sources/gaia/point_gaia.f
index 1781ba68b17a980d997a3f9dce55b1ed900a778a..864f2f6b903b8a64ae2ccaaebbbae7c57ec80d71 100644
--- a/sources/gaia/point_gaia.f
+++ b/sources/gaia/point_gaia.f
@@ -646,6 +646,10 @@
 19    FORMAT(1X,'ERROR IN VARSOR: CHECK_NVAR IS',1X,I2,
      &     /,1X,'WHILE IT SHOULD BE EQUAL TO',1X,I2,
      &     /,'CHECK THE CORRESPONDENCE BETWEEN TEXTE/MNEMO AND VARSOR')
+      ALLOCATE(ISDIR_GAIA(MAXVAR))
+      ISDIR_GAIA = .FALSE.
+      !WAVE DIRECTION NEED SPECIAL TREATEMENT
+      ISDIR_GAIA(13)  = .TRUE.
 !
 !-----------------------------------------------------------------------
 ! !JAJ #### IF REQUIRED, HERE WE CAN READ THE INPUT SECTIONS FILE
diff --git a/sources/tomawac/deall_tomawac.f b/sources/tomawac/deall_tomawac.f
index e041737a248f90ef2db87299c707ee4e238c6495..8d2549fad1299983fe3c8883fca8034f66408b60 100644
--- a/sources/tomawac/deall_tomawac.f
+++ b/sources/tomawac/deall_tomawac.f
@@ -385,6 +385,7 @@
 !     OUT_HISTORY
 !
       CALL OUTHIST_END()
+      DEALLOCATE(ISDIR_TOM)
 !
 !***********************************************************************
 !
diff --git a/sources/tomawac/declarations_tomawac.f b/sources/tomawac/declarations_tomawac.f
index b8d1f4ecb47fa0a4f802327351cd5f96dfc4e1c2..f4c210a1f223039f5b18ed37104ea0f63a96d2ad 100644
--- a/sources/tomawac/declarations_tomawac.f
+++ b/sources/tomawac/declarations_tomawac.f
@@ -953,6 +953,10 @@
 !
       LOGICAL :: CALC_DISSIP
 !
+!     WHETHER AN OUTPUT VARIABLE IS A DIRECTION
+!
+      LOGICAL, ALLOCATABLE, TARGET :: ISDIR_TOM(:)
+!
 !     TITLE
 !
       CHARACTER (LEN=80) :: TITCAS
diff --git a/sources/tomawac/point_tomawac.f b/sources/tomawac/point_tomawac.f
index 75c37696e99f0d44b6630a6a2ca2a8fa1c859005..bf04bd9c0ecab691501e2114d9fd6afa1c375e49 100644
--- a/sources/tomawac/point_tomawac.f
+++ b/sources/tomawac/point_tomawac.f
@@ -603,6 +603,11 @@
       CALL ADDBLO(VARSOR,DISSIP_ROL)
 !     40: PEAK DIRECTION
       CALL ADDBLO(VARSOR,SDPIC)
+!
+      ALLOCATE(ISDIR_TOM(MAXVAR))
+      ISDIR_TOM     = .FALSE.
+      ISDIR_TOM(3)  = .TRUE.
+      ISDIR_TOM(40) = .TRUE.
 
 !
 !.....BLOCK FOR VALIDATION
diff --git a/sources/tomawac/wac.F b/sources/tomawac/wac.F
index cc5e23b65f0d1551ed9713bf7351158b18d7379c..c90c04ed6466d47d394be68d8f9bd0f9f39cd0d8 100644
--- a/sources/tomawac/wac.F
+++ b/sources/tomawac/wac.F
@@ -494,7 +494,7 @@
 !         TIME SERIES FILE
           IF(YAHIST) THEN
             CALL OUTHIST('WAC   ',WAC_FILES(WACHI2),AT,NPOIN2,1,MAXVAR,
-     &                   SORLEO,VARSOR)
+     &                   SORLEO,VARSOR,ISDIR_TOM)
           ENDIF
           IF(IMPRES) THEN
 !
diff --git a/sources/tomawac/wac_init.F b/sources/tomawac/wac_init.F
index 496f6c5383b713721ae0f0011f777b369deb8e23..313ce0dfcc6d533283dcf38ec270a9f831242903 100644
--- a/sources/tomawac/wac_init.F
+++ b/sources/tomawac/wac_init.F
@@ -492,7 +492,7 @@
 !
       IF(YAHIST) THEN
         CALL OUTHIST('WAC   ',WAC_FILES(WACHI2),AT,NPOIN2,1,MAXVAR,
-     &               SORLEO,VARSOR)
+     &               SORLEO,VARSOR,ISDIR_TOM)
       ENDIF
 ! ---------------------------------------------------------------------
 !
diff --git a/sources/utils/bief/out_history.F90 b/sources/utils/bief/out_history.F90
index 4bb23ed93d389abab5854cc23540c6bc233a00da..b9d431f90096ab92df8e3aa0cf7d582808ad8089 100644
--- a/sources/utils/bief/out_history.F90
+++ b/sources/utils/bief/out_history.F90
@@ -743,7 +743,7 @@
                     SUBROUTINE OUTHIST &
 !                   ******************
 !
-      (MODNAME,OUTFILE,TIME,NPOIN2,NPLAN,MAXVAR,SORGND,VARSOR)
+      (MODNAME,OUTFILE,TIME,NPOIN2,NPLAN,MAXVAR,SORGND,VARSOR,ISDIR)
 !
 !***********************************************************************
 ! BIEF
@@ -769,6 +769,7 @@
 !>@ param[in]   NPLAN     Number of vertical planes
 !>@ param[in]   MAXVAR    Maximum number of variables
 !>@ param[in]   SORGND    Mask specifying which variable to write
+!>@ param[in]   ISDIR     Mask specifying which variable is a direction
 !>@ param[in]   VARSOR    Data of the variables for output
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !
@@ -781,13 +782,14 @@
       INTEGER, INTENT(IN) :: MAXVAR
       INTEGER, INTENT(IN) :: NPOIN2,NPLAN
       LOGICAL, INTENT(IN) :: SORGND(MAXVAR)
+      LOGICAL, INTENT(IN), OPTIONAL :: ISDIR(MAXVAR)
       TYPE(BIEF_OBJ), INTENT(IN) :: VARSOR
 !
 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
 !
       INTEGER IFILE,I,IERR,IPLAN,IVAR,NRCELLLOC,NRPOINT,HISTSTEP
       DOUBLE PRECISION, POINTER, DIMENSION(:) :: TMPDAT, OUTDAT
-      LOGICAL LEO
+      LOGICAL LEO,INTDIR
 !
 !-----------------------------------------------------------------------
 !
@@ -822,6 +824,11 @@
       IVAR = 0
       DO I=1,MAXVAR
         IF(SORGND(I)) THEN
+          IF(PRESENT(ISDIR)) THEN
+            INTDIR = ISDIR(I)
+          ELSE
+            INTDIR = .FALSE.
+          ENDIF
           IVAR = IVAR +1
           LOCDAT(IFILE)%DATAARR = 0.D0
           ! Interpolation of 2-d variables
@@ -829,7 +836,7 @@
             TMPDAT => VARSOR%ADR(I)%P%R((IPLAN-1)*NPOIN2+1:IPLAN*NPOIN2)
             OUTDAT => LOCDAT(IFILE)%DATAARR((IPLAN-1)*NRCELLLOC+1:IPLAN*NRCELLLOC)
             CALL INTERP_POINT(INTDAT(IFILE),NRPOINT,TMPDAT,OUTDAT, &
-                              NRCELLLOC)
+                              NRCELLLOC,INTDIR)
           ENDDO
           ! Add data to file
           CALL ADD_DATA(OUTFILE%FMT,OUTFILE%LU, &
@@ -1190,7 +1197,7 @@
                     SUBROUTINE INTERP_POINT &
 !                   ***********************
 !
-      (INTVAR,NRPOINTS,DATA_IN,DATA_OUT,OUT_SIZE)
+      (INTVAR,NRPOINTS,DATA_IN,DATA_OUT,OUT_SIZE,ISDIR)
 !
 !***********************************************************************
 ! BIEF
@@ -1210,6 +1217,10 @@
 !>@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
+!>@param[in]      ISDIR    Indicates whether data contains directions
+!                          Directions are assumed to be in degrees
+!                          The method works for nautical and cartesian convention
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !
 !
@@ -1217,6 +1228,7 @@
 !
       TYPE(INTERPDATA), INTENT(IN) :: INTVAR
       INTEGER, INTENT(IN) :: NRPOINTS,OUT_SIZE
+      LOGICAL, INTENT(IN) :: ISDIR
       DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: DATA_IN
       DOUBLE PRECISION, DIMENSION(OUT_SIZE), INTENT(INOUT):: DATA_OUT
 !
@@ -1224,6 +1236,10 @@
 !
       INTEGER I,I1,I2,I3
       DOUBLE PRECISION, DIMENSION(NRPOINTS) :: TEMP_ARR
+      DOUBLE PRECISION COSDIR,SINDIR
+      DOUBLE PRECISION, PARAMETER :: PI = 4.D0 * ATAN( 1.D0 )
+      DOUBLE PRECISION, PARAMETER :: RADDEG = 180.D0/PI
+      DOUBLE PRECISION, PARAMETER :: DEGRAD = PI/180.D0
 !
 !-----------------------------------------------------------------------
 !
@@ -1238,9 +1254,19 @@
         ! Check for valid point
         IF(I1.GT.0 .AND. I2.GT.0 .AND. I3.GT.0) THEN
         ! Perform interpolation by matrix vector product
-          TEMP_ARR(I) = INTVAR%WEIGHTS(I,1)*DATA_IN(I1) +  &
-                        INTVAR%WEIGHTS(I,2)*DATA_IN(I2) +  &
-                        INTVAR%WEIGHTS(I,3)*DATA_IN(I3)
+          IF (.NOT.ISDIR) THEN
+            TEMP_ARR(I) = INTVAR%WEIGHTS(I,1)*DATA_IN(I1) +  &
+                          INTVAR%WEIGHTS(I,2)*DATA_IN(I2) +  &
+                          INTVAR%WEIGHTS(I,3)*DATA_IN(I3)
+          ELSE
+            COSDIR = INTVAR%WEIGHTS(I,1)*COS(DEGRAD*DATA_IN(I1)) +  &
+                     INTVAR%WEIGHTS(I,2)*COS(DEGRAD*DATA_IN(I2)) +  &
+                     INTVAR%WEIGHTS(I,3)*COS(DEGRAD*DATA_IN(I3))
+            SINDIR = INTVAR%WEIGHTS(I,1)*SIN(DEGRAD*DATA_IN(I1)) +  &
+                     INTVAR%WEIGHTS(I,2)*SIN(DEGRAD*DATA_IN(I2)) +  &
+                     INTVAR%WEIGHTS(I,3)*SIN(DEGRAD*DATA_IN(I3))
+            TEMP_ARR(I) = RADDEG*ATAN2(SINDIR,COSDIR)
+          ENDIF
         ENDIF
       ENDDO
       ! Copy data to right array