From 92c5994fffa215690d090669d89c92a4a6e7b579 Mon Sep 17 00:00:00 2001
From: Thomas SAILLOUR <saillour.thomas@gmail.com>
Date: Tue, 28 Jan 2025 13:40:10 +0000
Subject: [PATCH] [fix] Time Series extraction for longitude/latitude meshes
across IDL (International Date Line)
---
NEWS.txt | 2 +
sources/tomawac/wac_init.F | 3 +-
sources/utils/bief/out_history.F90 | 82 +++++++++++++++++++++++++-----
3 files changed, 73 insertions(+), 14 deletions(-)
diff --git a/NEWS.txt b/NEWS.txt
index 36126af6ff..eed12d2bb1 100644
--- a/NEWS.txt
+++ b/NEWS.txt
@@ -1,6 +1,8 @@
Latest changes
==============
+TELEMAC-2D/TELEMAC-3D/TOMAWAC/GAIA: fix time series exports for global models
+
MASCARET: New keyword VARIABLES PRECISION
to control the precision of saved variables
diff --git a/sources/tomawac/wac_init.F b/sources/tomawac/wac_init.F
index 44ba4ea807..c9f8d04a18 100644
--- a/sources/tomawac/wac_init.F
+++ b/sources/tomawac/wac_init.F
@@ -470,7 +470,8 @@
CALL OUTHIST_PREPARE('WAC ',WAC_FILES(WACCOO),
& WAC_FILES(WACHI2),MESH,NELEM2,NPOIN2,
- & DATE,TIME,MAXVAR,SORLEO,TEXTE,VARSOR,1)
+ & DATE,TIME,MAXVAR,SORLEO,TEXTE,VARSOR,1,
+ & SPHE=SPHE)
CALL HISTPERIOD('WAC ',AT,IMPRES_HIS)
IF(IMPRES.OR.IMPRES_HIS.OR.LT.EQ.0) THEN
diff --git a/sources/utils/bief/out_history.F90 b/sources/utils/bief/out_history.F90
index b069b3c120..b3bd98b8f1 100644
--- a/sources/utils/bief/out_history.F90
+++ b/sources/utils/bief/out_history.F90
@@ -365,7 +365,7 @@
! **************************
!
(MODNAME,COORFILE,OUTFILE,MESH2D,NELEM2,NPOIN2,DATE,TIME,MAXVAR, &
- SORG2D,TEXTE,VARSOR,NPLAN,PROTYP,LATIT,LONGIT,KEEP_LONLAT)
+ SORG2D,TEXTE,VARSOR,NPLAN,PROTYP,LATIT,LONGIT,KEEP_LONLAT,SPHE)
!
!***********************************************************************
! BIEF
@@ -403,6 +403,8 @@
!>@param[in] LONGIT Longitude of origin point (Mercator conversion)
!>@param[in] KEEP_LONLAT If YES and PROTYP = 3, results are written in
! longitude/latitude, not Mercator projection
+!>@param[in] SPHE If YES, TOMAWAC results are computed
+! in longitude/latitude (SPHERICAL COORDINATES)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
USE DECLARATIONS_PARALLEL
@@ -421,7 +423,7 @@
TYPE(BIEF_OBJ), INTENT(IN) :: VARSOR
INTEGER, OPTIONAL, INTENT(IN) :: PROTYP
DOUBLE PRECISION, OPTIONAL, INTENT(IN) :: LATIT,LONGIT
- LOGICAL, OPTIONAL, INTENT(IN) :: KEEP_LONLAT
+ LOGICAL, OPTIONAL, INTENT(IN) :: KEEP_LONLAT,SPHE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
@@ -432,7 +434,7 @@
DOUBLE PRECISION, DIMENSION(3) :: DATADD3
DOUBLE PRECISION, PARAMETER :: R=6.37D6
DOUBLE PRECISION PS4,TAN1,TAN2,LONGIRAD,LATI0,LONGI0
- LOGICAL YESLL
+ LOGICAL YESLL,SPHERI
! Interpolation method: set to 1 for linear interpolation or 0 for
! nearest neighbour interpolation
@@ -468,6 +470,12 @@
ELSE
YESLL = .TRUE.
ENDIF
+!
+ IF(PRESENT(SPHE)) THEN
+ SPHERI = SPHE
+ ELSE
+ SPHERI = .FALSE.
+ ENDIF
!
!-----------------------------------------------------------------------
!
@@ -635,7 +643,7 @@
CALL INTERP_POINT_PREPARE(TCDAT(IFILE)%XHISTOUT, &
TCDAT(IFILE)%YHISTOUT,TCDAT(IFILE)%NRPOINT,INTDAT(IFILE), &
- INTERP_METHOD,MESH2D)
+ INTERP_METHOD,MESH2D,MODNAME,PROJECTION,SPHERI)
! Determine which points to consider
INTDAT(IFILE)%MASKARR = INTDAT(IFILE)%POINTLIST(:,1).NE.0
@@ -1042,7 +1050,7 @@
SUBROUTINE INTERP_POINT_PREPARE &
! *******************************
!
- (XINT,YINT,NPINT,INTVAR,INTERP_METHOD,MESH2D)
+ (XINT,YINT,NPINT,INTVAR,INTERP_METHOD,MESH2D,MODNAME,PROTYP,SPHERI)
!
!***********************************************************************
! BIEF
@@ -1066,6 +1074,10 @@
!>@param[in] INTERP_METHOD 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
+!>@param[in] PROTYP Spatial projection type (3: longitude/latitude)
+!>@param[in] SPHERI If YES, TOMAWAC results are computed
+! in longitude/latitude (SPHERICAL COORDINATES)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!
@@ -1077,6 +1089,9 @@
INTEGER, INTENT(IN) :: INTERP_METHOD
TYPE(BIEF_MESH), INTENT(IN) :: MESH2D
+ CHARACTER(LEN=6), INTENT(IN) :: MODNAME
+ INTEGER, INTENT(IN) :: PROTYP
+ LOGICAL, INTENT(IN) :: SPHERI
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
@@ -1116,7 +1131,7 @@
YTRI(3) = MESH2D%Y%R(IKLE2(J3))
! Check if a point is in a trinagle
- CALL IN_TRIANGLE(XP,YP,XTRI,YTRI,INTRI)
+ CALL IN_TRIANGLE(XP,YP,XTRI,YTRI,INTRI,MODNAME,PROTYP,SPHERI)
IF(INTRI) THEN
! Store points
@@ -1239,7 +1254,7 @@
SUBROUTINE IN_TRIANGLE &
! **********************
!
- (XP,YP,XTRI,YTRI,INTRI)
+ (XP,YP,XTRI,YTRI,INTRI,MODNAME,PROTYP,SPHERI)
!
!***********************************************************************
! BIEF
@@ -1260,6 +1275,10 @@
!>@param[in] XTRI x-coordinates of the triangle
!>@param[in] YTRI y-coordinates of the triangle
!>@param[in,out] INTRI true if the point is in a triangle
+!>@param [in] MODNAME Name of the calling module
+!>@param[in] PROTYP Spatial projection type (3: longitude/latitude)
+!>@param[in] SPHERI If YES, TOMAWAC results are computed
+! in longitude/latitude (SPHERICAL COORDINATES)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!
@@ -1268,21 +1287,58 @@
DOUBLE PRECISION, INTENT(IN) :: XP,YP
DOUBLE PRECISION, INTENT(IN), DIMENSION(3) :: XTRI,YTRI
LOGICAL, INTENT(INOUT) :: INTRI
+ CHARACTER(LEN=6), INTENT(IN) :: MODNAME
+ INTEGER, INTENT(IN) :: PROTYP
+ LOGICAL, INTENT(IN) :: SPHERI
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
DOUBLE PRECISION, DIMENSION(2) :: V0,V1,V2
DOUBLE PRECISION DOT00,DOT01,DOT02,DOT11,DOT12,INV_DENOM,U,V
+ DOUBLE PRECISION R,PIR,TWOPIR
+ LOGICAL :: IDLCORRECT = .FALSE.
!
!-----------------------------------------------------------------------
!
+ IF(MODNAME(1:3).EQ.'WAC') THEN
+ PIR = 180.D0
+ ELSE
+ R = 6370000.D0
+ PIR = 4.D0 * ATAN(1.D0) * R
+ ENDIF
+ TWOPIR = 2.D0 * PIR
+
+ ! Correct the coordinates on either parts of the IDL
+ ! only for MERCATOR (2) or LAT/LON (3) projections
+ IF(PROTYP.GT.1.OR.(MODNAME(1:3).EQ.'WAC'.AND.SPHERI)) THEN
+ IDLCORRECT = .TRUE.
+ ENDIF
+
! Compute vectors
- V0(1) = XTRI(3) - XTRI(1)
- V0(2) = YTRI(3) - YTRI(1)
- V1(1) = XTRI(2) - XTRI(1)
- V1(2) = YTRI(2) - YTRI(1)
- V2(1) = XP - XTRI(1)
- V2(2) = YP - YTRI(1)
+ V0(1) = XTRI(3) - XTRI(1)
+ V0(2) = YTRI(3) - YTRI(1)
+ V1(1) = XTRI(2) - XTRI(1)
+ V1(2) = YTRI(2) - YTRI(1)
+ V2(1) = XP - XTRI(1)
+ V2(2) = YP - YTRI(1)
+
+ IF(IDLCORRECT) THEN
+ IF(V0(1).GT.PIR) THEN
+ V0(1) = V0(1) - TWOPIR
+ ELSEIF(V0(1).LT.-PIR) THEN
+ V0(1) = V0(1) + TWOPIR
+ ENDIF
+ IF(V1(1).GT.PIR) THEN
+ V1(1) = V1(1) - TWOPIR
+ ELSEIF(V1(1).LT.-PIR) THEN
+ V1(1) = V1(1) + TWOPIR
+ ENDIF
+ IF(V2(1).GT.PIR) THEN
+ V2(1) = V2(1) - TWOPIR
+ ELSEIF(V2(1).LT.-PIR) THEN
+ V2(1) = V2(1) + TWOPIR
+ ENDIF
+ ENDIF
! Compute dot products
DOT00 = DOT_PRODUCT(V0,V0)
--
GitLab