Philippe 08/04/2022: stations: add position and interpolation coefficients in TSTATIONDATA type

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'

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)

src/MNH/station_tools.f90

+!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

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