From ae50524e13a2b085ec398c859d38f07fa4cbe275 Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Fri, 8 Apr 2022 11:51:10 +0200
Subject: [PATCH] Philippe 08/04/2022: stations: add position and interpolation
coefficients in TSTATIONDATA type
---
src/MNH/ini_surfstationn.f90 | 3 ++
src/MNH/modd_type_station.f90 | 15 +++++-
src/MNH/station_tools.f90 | 95 +++++++++++++++++++++++++++++++++++
src/MNH/stationn.f90 | 43 ++++++++--------
4 files changed, 134 insertions(+), 22 deletions(-)
create mode 100644 src/MNH/station_tools.f90
diff --git a/src/MNH/ini_surfstationn.f90 b/src/MNH/ini_surfstationn.f90
index f68f35028..3cde34d16 100644
--- a/src/MNH/ini_surfstationn.f90
+++ b/src/MNH/ini_surfstationn.f90
@@ -262,6 +262,8 @@ END SUBROUTINE ALLOCATE_STATION_n
!----------------------------------------------------------------------------
SUBROUTINE INI_INTERP_STATION_n()
!
+USE MODE_STATION_TOOLS, ONLY: STATION_POSITION
+
INTEGER :: JII
INTEGER :: IIU, IJU
!
@@ -271,6 +273,7 @@ IF ( ALL(TSTATIONS(:)%XLAT /= XUNDEF) .AND. ALL(TSTATIONS(:)%XLON /= XUNDEF) ) T
CALL SM_XYHAT(PLATOR,PLONOR, &
TSTATIONS(JII)%XLAT, TSTATIONS(JII)%XLON, &
TSTATIONS(JII)%XX, TSTATIONS(JII)%XY )
+ CALL STATION_POSITION( TSTATIONS(JII) )
END DO
ELSE
CMNHMSG(1) = 'Error in station position'
diff --git a/src/MNH/modd_type_station.f90 b/src/MNH/modd_type_station.f90
index b929d28d4..ac0134eab 100644
--- a/src/MNH/modd_type_station.f90
+++ b/src/MNH/modd_type_station.f90
@@ -37,7 +37,7 @@
! ------------
!
use modd_type_date, only: date_time
-use modd_parameters, only: NUNDEF, XUNDEF
+use modd_parameters, only: NNEGUNDEF, NUNDEF, XUNDEF
implicit none
@@ -57,6 +57,7 @@ TYPE TSTATIONDATA
CHARACTER(LEN=8) :: CNAME = '' ! station name
CHARACTER(LEN=8) :: CTYPE = '' ! station type (currently not used)
LOGICAL :: LERROR = .FALSE. !
+LOGICAL :: LPRESENT = .FALSE. ! If true, this station is situated on this process
REAL :: XX = XUNDEF ! X(n)
REAL :: XY = XUNDEF ! Y(n)
@@ -65,6 +66,18 @@ REAL :: XLON = XUNDEF ! longitude(n)
REAL :: XLAT = XUNDEF ! latitude (n)
REAL :: XZS = XUNDEF ! zs(n)
+! Position in the mesh
+INTEGER :: NI_M = NNEGUNDEF ! X position for mass-point axis (between this one and the next one)
+INTEGER :: NJ_M = NNEGUNDEF ! Y position for mass-point axis (between this one and the next one)
+INTEGER :: NI_U = NNEGUNDEF ! X position for u-point axis (between this one and the next one)
+INTEGER :: NJ_V = NNEGUNDEF ! Y position for v-point axis (between this one and the next one)
+
+! Coefficient to interpolate values (stations are usually not exactly on mesh points)
+REAL :: XXMCOEF = XUNDEF ! Interpolation coefficient for X (mass-point)
+REAL :: XYMCOEF = XUNDEF ! Interpolation coefficient for Y (mass-point)
+REAL :: XXUCOEF = XUNDEF ! Interpolation coefficient for X (U-point)
+REAL :: XYVCOEF = XUNDEF ! Interpolation coefficient for Y (V-point)
+
INTEGER :: NK = NUNDEF ! Model level for altitude comparisons
REAL, DIMENSION(:), ALLOCATABLE :: XZON ! zonal wind(n)
diff --git a/src/MNH/station_tools.f90 b/src/MNH/station_tools.f90
new file mode 100644
index 000000000..1aa02f68e
--- /dev/null
+++ b/src/MNH/station_tools.f90
@@ -0,0 +1,95 @@
+!MNH_LIC Copyright 2022-2022 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+! Author:
+! P. Wautelet 08/04/2022
+!-----------------------------------------------------------------
+! Modifications:
+!-----------------------------------------------------------------
+! ##################
+MODULE MODE_STATION_TOOLS
+! ##################
+
+IMPLICIT NONE
+
+PRIVATE
+
+PUBLIC :: STATION_POSITION
+
+CONTAINS
+
+! ######################################
+SUBROUTINE STATION_POSITION( TPSTATION )
+! ######################################
+! Subroutine to determine the position of a station on the model grid
+! and set the useful coefficient for data interpolation
+
+ USE MODD_CONF, ONLY: L1D
+ USE MODD_GRID_n, ONLY: XXHAT, XYHAT
+ USE MODD_TYPE_STATION, ONLY: TSTATIONDATA
+
+ USE MODE_TOOLS_ll, ONLY: GET_INDICE_ll, LWEST_ll, LEAST_ll, LNORTH_ll, LSOUTH_ll
+
+ IMPLICIT NONE
+
+ TYPE(TSTATIONDATA), INTENT(INOUT) :: TPSTATION
+
+ INTEGER :: IIB ! domain sizes of current process
+ INTEGER :: IJB !
+ INTEGER :: IIE !
+ INTEGER :: IJE !
+ INTEGER :: IIU !
+ INTEGER :: IJU !
+ REAL, DIMENSION(SIZE(XXHAT)) :: ZXHATM ! mass point coordinates
+ REAL, DIMENSION(SIZE(XYHAT)) :: ZYHATM ! mass point coordinates
+
+ IIU = SIZE( XXHAT )
+ IJU = SIZE( XYHAT )
+
+ CALL GET_INDICE_ll (IIB, IJB, IIE, IJE )
+
+ ! Interpolations of model variables to mass points
+ ZXHATM(1:IIU-1) = 0.5 * XXHAT(1:IIU-1) + 0.5 * XXHAT(2:IIU )
+ ZXHATM( IIU ) = 1.5 * XXHAT( IIU ) - 0.5 * XXHAT( IIU-1)
+
+ ZYHATM(1:IJU-1) = 0.5 * XYHAT(1:IJU-1) + 0.5 * XYHAT(2:IJU )
+ ZYHATM( IJU ) = 1.5 * XYHAT( IJU ) - 0.5 * XYHAT( IJU-1)
+
+ TPSTATION%LPRESENT = .FALSE.
+
+ ! X position
+ TPSTATION%NI_U = COUNT( XXHAT (:) <= TPSTATION%XX )
+ TPSTATION%NI_M = COUNT( ZXHATM(:) <= TPSTATION%XX )
+
+ IF ( TPSTATION%NI_M<=IIB-1 .AND. LWEST_ll() .AND. .NOT. L1D ) TPSTATION%LERROR = .TRUE.
+ IF ( TPSTATION%NI_M>=IIE .AND. LEAST_ll() .AND. .NOT. L1D ) TPSTATION%LERROR = .TRUE.
+
+ ! Y position
+ TPSTATION%NJ_V = COUNT( XYHAT (:) <= TPSTATION%XY )
+ TPSTATION%NJ_M = COUNT( ZYHATM(:) <= TPSTATION%XY )
+
+ IF ( TPSTATION%NJ_M<=IJB-1 .AND. LSOUTH_ll() .AND. .NOT. L1D ) TPSTATION%LERROR = .TRUE.
+ IF ( TPSTATION%NJ_M>=IJE .AND. LNORTH_ll() .AND. .NOT. L1D ) TPSTATION%LERROR = .TRUE.
+
+ ! Position of station according to processes
+ IF ( TPSTATION%NI_U >= IIB .AND. TPSTATION%NI_U <= IIE &
+ .AND. TPSTATION%NJ_V >= IJB .AND. TPSTATION%NJ_V <= IJE ) TPSTATION%LPRESENT = .TRUE.
+ IF ( L1D ) TPSTATION%LPRESENT = .TRUE.
+
+ ! Computations only on correct process
+ IF ( TPSTATION%LPRESENT .AND. .NOT. L1D ) THEN
+ ! Interpolation coefficient for X (mass-point)
+ TPSTATION%XXMCOEF = ( TPSTATION%XX - ZXHATM(TPSTATION%NI_M) ) / ( ZXHATM(TPSTATION%NI_M+1) - ZXHATM(TPSTATION%NI_M) )
+ ! Interpolation coefficient for Y (mass-point)
+ TPSTATION%XYMCOEF = ( TPSTATION%XY - ZYHATM(TPSTATION%NJ_M) ) / ( ZYHATM(TPSTATION%NJ_M+1) - ZYHATM(TPSTATION%NJ_M) )
+ ! Interpolation coefficient for X (U-point)
+ TPSTATION%XXUCOEF = ( TPSTATION%XX - XXHAT(TPSTATION%NI_U) ) / ( XXHAT(TPSTATION%NI_U+1) - XXHAT(TPSTATION%NI_U) )
+ ! Interpolation coefficient for Y (V-point)
+ TPSTATION%XYVCOEF = ( TPSTATION%XY - XYHAT(TPSTATION%NJ_V) ) / ( XYHAT(TPSTATION%NJ_V+1) - XYHAT(TPSTATION%NJ_V) )
+ END IF
+
+END SUBROUTINE STATION_POSITION
+
+END MODULE MODE_STATION_TOOLS
diff --git a/src/MNH/stationn.f90 b/src/MNH/stationn.f90
index 92105242e..d59477e66 100644
--- a/src/MNH/stationn.f90
+++ b/src/MNH/stationn.f90
@@ -152,7 +152,6 @@ REAL :: ZU_STAT ! horizontal wind speed at station location (along x)
REAL :: ZV_STAT ! horizontal wind speed at station location (along y)
REAL :: ZGAM ! rotation between meso-nh base and spherical lat-lon base.
!
-INTEGER :: IINFO_ll ! return code
INTEGER :: IRESP ! return code
INTEGER :: I ! loop for stations
INTEGER :: J ! loop for levels
@@ -318,7 +317,7 @@ END IF
IF (GSTORE) THEN
DO I=1,NUMBSTAT
!
- IF ((ZTHIS_PROCS(I)==1.).AND.(.NOT. TSTATIONS(I)%LERROR)) THEN
+ IF ( TSTATIONS(I)%LPRESENT .AND. .NOT. TSTATIONS(I)%LERROR ) THEN
IF (TSTATIONS(I)%NK/= XUNDEF) THEN
J = TSTATIONS(I)%NK
ELSE ! suppose TSTATIONS(I)%XZ /= XUNDEF
@@ -461,17 +460,17 @@ ELSEIF (L1D) THEN
JI=2
JJ=2
ELSE
- JI=II(I)
- JJ=IJ(I)
+ JI=TSTATIONS(I)%NI_M
+ JJ=TSTATIONS(I)%NJ_M
END IF
!
!
IF ((JI .GE. 1).AND. (JI .LE. SIZE(PA,1)) .AND. &
(JJ .GE. 1).AND. (JJ .LE. SIZE(PA,2))) &
-PB = (1.-ZYCOEF(I)) * (1.-ZXCOEF(I)) * PA(JI,JJ) + &
- (1.-ZYCOEF(I)) * (ZXCOEF(I)) * PA(JI+1,JJ) + &
- ( ZYCOEF(I)) * (1.-ZXCOEF(I)) * PA(JI,JJ+1) + &
- ( ZYCOEF(I)) * (ZXCOEF(I)) * PA(JI+1,JJ+1)
+PB = (1.-TSTATIONS(I)%XYMCOEF) * (1.-TSTATIONS(I)%XXMCOEF) * PA(JI,JJ) + &
+ (1.-TSTATIONS(I)%XYMCOEF) * (TSTATIONS(I)%XXMCOEF) * PA(JI+1,JJ) + &
+ ( TSTATIONS(I)%XYMCOEF) * (1.-TSTATIONS(I)%XXMCOEF) * PA(JI,JJ+1) + &
+ ( TSTATIONS(I)%XYMCOEF) * (TSTATIONS(I)%XXMCOEF) * PA(JI+1,JJ+1)
!
END FUNCTION STATION_INTERP_2D
!----------------------------------------------------------------------------
@@ -491,16 +490,16 @@ ELSEIF (L1D) THEN
JI=2
JJ=2
ELSE
- JI=II(I)
- JJ=IJ(I)
+ JI=TSTATIONS(I)%NI_M
+ JJ=TSTATIONS(I)%NJ_M
END IF
!
IF ((JI .GE. 1).AND. (JI .LE. SIZE(PA,1)) .AND. &
(JJ .GE. 1).AND. (JJ .LE. SIZE(PA,2))) &
-PB = (1.- ZYCOEF(I)) * (1.-ZUCOEF(I)) * PA(JI ,JJ ) &
- + (1.- ZYCOEF(I)) * ( ZUCOEF(I)) * PA(JI+1,JJ ) &
- + ( ZYCOEF(I)) * (1.-ZUCOEF(I)) * PA(JI ,JJ+1) &
- + ( ZYCOEF(I)) * ( ZUCOEF(I)) * PA(JI+1,JJ+1)
+PB = (1.- TSTATIONS(I)%XYMCOEF) * (1.-TSTATIONS(I)%XXUCOEF) * PA(JI ,JJ ) &
+ + (1.- TSTATIONS(I)%XYMCOEF) * ( TSTATIONS(I)%XXUCOEF) * PA(JI+1,JJ ) &
+ + ( TSTATIONS(I)%XYMCOEF) * (1.-TSTATIONS(I)%XXUCOEF) * PA(JI ,JJ+1) &
+ + ( TSTATIONS(I)%XYMCOEF) * ( TSTATIONS(I)%XXUCOEF) * PA(JI+1,JJ+1)
!
END FUNCTION STATION_INTERP_2D_U
!----------------------------------------------------------------------------
@@ -520,16 +519,16 @@ ELSEIF (L1D) THEN
JI=2
JJ=2
ELSE
- JI=II(I)
- JJ=IJ(I)
+ JI=TSTATIONS(I)%NI_M
+ JJ=TSTATIONS(I)%NJ_M
END IF
!
IF ((JI .GT. 0).AND. (JI .LT. SIZE(PA,1)) .AND. &
(JJ .GT. 0).AND. (JJ .LT. SIZE(PA,2))) &
-PB = (1.- ZVCOEF(I)) * (1.-ZXCOEF(I)) * PA(JI ,JJ ) &
- + (1.- ZVCOEF(I)) * ( ZXCOEF(I)) * PA(JI+1,JJ ) &
- + ( ZVCOEF(I)) * (1.-ZXCOEF(I)) * PA(JI ,JJ+1) &
- + ( ZVCOEF(I)) * ( ZXCOEF(I)) * PA(JI+1,JJ+1)
+PB = (1.- TSTATIONS(I)%XYVCOEF) * (1.-TSTATIONS(I)%XXMCOEF) * PA(JI ,JJ ) &
+ + (1.- TSTATIONS(I)%XYVCOEF) * ( TSTATIONS(I)%XXMCOEF) * PA(JI+1,JJ ) &
+ + ( TSTATIONS(I)%XYVCOEF) * (1.-TSTATIONS(I)%XXMCOEF) * PA(JI ,JJ+1) &
+ + ( TSTATIONS(I)%XYVCOEF) * ( TSTATIONS(I)%XXMCOEF) * PA(JI+1,JJ+1)
!
END FUNCTION STATION_INTERP_2D_V
!----------------------------------------------------------------------------
@@ -538,7 +537,9 @@ SUBROUTINE DISTRIBUTE_STATION(PAS)
!
REAL, INTENT(INOUT) :: PAS
!
-PAS = PAS * ZTHIS_PROCS(I)
+INTEGER :: IINFO_ll ! return code
+
+IF ( .NOT. TSTATIONS(I)%LPRESENT ) PAS = 0.
CALL REDUCESUM_ll(PAS,IINFO_ll)
!
END SUBROUTINE DISTRIBUTE_STATION
--
GitLab