From 55f55971dcf5c56ccedf6a1cc6b44665657faffc Mon Sep 17 00:00:00 2001
From: Alexander BREUGEM <alexander.breugem@imdc.be>
Date: Thu, 6 Feb 2025 10:10:33 +0000
Subject: [PATCH] Fix the interpolation of wave directions in time series files
---
NEWS.txt | 2 +
.../tomawac/3Dcoupling/coordinate_tom.txt | 22 ++++++++++
examples/tomawac/3Dcoupling/fom_histcoup.slf | 3 ++
.../tomawac/3Dcoupling/tom_3Dcouplittoral.cas | 2 +
.../tomawac/3Dcoupling/vnv_3d_coupling.py | 12 ++++++
sources/gaia/deall_gaia.f | 1 +
sources/gaia/declarations_gaia.f | 4 ++
sources/gaia/gaia_init.F | 2 +-
sources/gaia/gaia_write_results.f | 2 +-
sources/gaia/point_gaia.f | 4 ++
sources/tomawac/deall_tomawac.f | 1 +
sources/tomawac/declarations_tomawac.f | 4 ++
sources/tomawac/point_tomawac.f | 5 +++
sources/tomawac/wac.F | 2 +-
sources/tomawac/wac_init.F | 2 +-
sources/utils/bief/out_history.F90 | 40 +++++++++++++++----
16 files changed, 97 insertions(+), 11 deletions(-)
create mode 100644 examples/tomawac/3Dcoupling/coordinate_tom.txt
create mode 100644 examples/tomawac/3Dcoupling/fom_histcoup.slf
diff --git a/NEWS.txt b/NEWS.txt
index faeab047ff..8e553e71ae 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 0000000000..cfa3d5d87d
--- /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 0000000000..fc4173b0a2
--- /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 93f9102adc..979e8adff8 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 167c23f5b8..94e9ba7dbe 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 de02118806..7a8c30fe8d 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 9e187902ab..35eb291b1c 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 9e89f6da39..9d5678c29d 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 01fcaba5db..4f6582a528 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 1781ba68b1..864f2f6b90 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 e041737a24..8d2549fad1 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 b8d1f4ecb4..f4c210a1f2 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 75c37696e9..bf04bd9c0e 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 cc5e23b65f..c90c04ed64 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 496f6c5383..313ce0dfcc 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 4bb23ed93d..b9d431f900 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
--
GitLab