From 488baaf1adbb7f14b8298971463f13809f3e517e Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Wed, 20 Apr 2022 11:35:49 +0200
Subject: [PATCH] Philippe 20/04/2022: stations: optimise stations (less
communications, less memory...) Interpolation was not correct for U and V
(did not use correct position indices)
---
src/MNH/default_desfmn.f90 | 3 +-
src/MNH/ini_modeln.f90 | 6 +-
src/MNH/ini_surfstationn.f90 | 295 +++++++++-----------------
src/MNH/modd_allstationn.f90 | 16 +-
src/MNH/modd_stationn.f90 | 18 +-
src/MNH/modd_type_station.f90 | 12 +-
src/MNH/modn_stationn.f90 | 10 +-
src/MNH/station_reader.f90 | 162 ++++++++-------
src/MNH/station_tools.f90 | 376 ++++++++++++++++++++++++++++++----
src/MNH/stationn.f90 | 351 ++++++-------------------------
src/MNH/write_stationn.f90 | 228 +++++++++++++++++++--
11 files changed, 819 insertions(+), 658 deletions(-)
diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index c1be2c51e..c9f530786 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-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.
@@ -603,7 +603,6 @@ XZ_STAT(:) = XUNDEF
XLAT_STAT(:) = XUNDEF
XLON_STAT(:) = XUNDEF
CNAME_STAT(:) = ''
-CTYPE_STAT(:) = ''
CFILE_STAT = 'NO_INPUT_CSV'
LDIAG_SURFRAD = .TRUE.
!
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index 718f11ce6..6ee23d80a 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-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.
@@ -2565,9 +2565,7 @@ CALL INI_AIRCRAFT_BALLOON(TPINIFILE,XTSTEP, TDTSEG, XSEGLEN, NRR, NSV, &
!* 24. STATION initializations
! -----------------------
!
-CALL INI_SURFSTATION_n(XTSTEP, XSEGLEN, NRR, NSV, &
- CTURB=="TKEL" , KMI, &
- XLATORI, XLONORI )
+CALL INI_SURFSTATION_n( )
!
!-------------------------------------------------------------------------------
!
diff --git a/src/MNH/ini_surfstationn.f90 b/src/MNH/ini_surfstationn.f90
index 3cde34d16..ebd004817 100644
--- a/src/MNH/ini_surfstationn.f90
+++ b/src/MNH/ini_surfstationn.f90
@@ -9,21 +9,7 @@ MODULE MODI_INI_SURFSTATION_n
!
INTERFACE
!
- SUBROUTINE INI_SURFSTATION_n(PTSTEP, PSEGLEN, &
- KRR, KSV, OUSETKE, KMI, &
- PLATOR, PLONOR )
-!
-USE MODD_TYPE_DATE
-REAL, INTENT(IN) :: PTSTEP ! time step
-REAL, INTENT(IN) :: PSEGLEN ! segment length
-INTEGER, INTENT(IN) :: KRR ! number of moist variables
-INTEGER, INTENT(IN) :: KSV ! number of scalar variables
-LOGICAL, INTENT(IN) :: OUSETKE ! flag to use tke
-REAL, INTENT(IN) :: PLATOR ! latitude of origine point
-REAL, INTENT(IN) :: PLONOR ! longitude of origine point
-INTEGER, INTENT(IN) :: KMI ! MODEL NUMBER
-!
-!-------------------------------------------------------------------------------
+ SUBROUTINE INI_SURFSTATION_n( )
!
END SUBROUTINE INI_SURFSTATION_n
!
@@ -31,11 +17,9 @@ END INTERFACE
!
END MODULE MODI_INI_SURFSTATION_n
!
-! ########################################################
- SUBROUTINE INI_SURFSTATION_n(PTSTEP, PSEGLEN, &
- KRR, KSV, OUSETKE, KMI, &
- PLATOR, PLONOR )
-! ########################################################
+! ###############################
+ SUBROUTINE INI_SURFSTATION_n( )
+! ###############################
!
!
!!**** *INI_SURFSTATION_n* -
@@ -69,222 +53,143 @@ END MODULE MODI_INI_SURFSTATION_n
! P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
! R. Schoetter 11/2019: work for cartesian coordinates + parallel.
! E.Jezequel 02/2021: read stations from CVS file
-! P. Wautelet 07/04/2022: rewrite types for stations
+! P. Wautelet 04/2022: restructure stations for better performance, reduce memory usage and correct some problems/bugs
! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE MODD_ALLSTATION_n
-USE MODD_CONF
-USE MODD_DIM_n
-USE MODD_DYN_n
-USE MODD_GRID
-USE MODD_GRID_n
-USE MODD_NESTING
-USE MODD_PARAMETERS
-USE MODD_SHADOWS_n
+USE MODD_CONF, ONLY: LCARTESIAN
+USE MODD_DYN, ONLY: XSEGLEN
+USE MODD_DYN_n, ONLY: DYN_MODEL, XTSTEP
+USE MODD_GRID_n, ONLY: XXHAT, XYHAT
+USE MODD_PARAMETERS, ONLY: JPHEXT, JPVEXT
USE MODD_STATION_n
-USE MODD_TYPE_DATE
-USE MODD_VAR_ll, ONLY: IP
+USE MODD_TYPE_STATION
!
-USE MODE_GATHER_ll
-USE MODE_GRIDPROJ
-USE MODE_ll
USE MODE_MSG
-!
-USE MODI_INI_STATION_N
+USE MODE_STATION_READER
+USE MODE_STATION_TOOLS, ONLY: STATION_ADD, STATION_ALLOCATE, STATION_INI_INTERP, STATION_POSITION
+USE MODE_ALLOCBUFFER_ll, ONLY: ALLOCBUFFER_ll
+USE MODE_GATHER_ll, ONLY: GATHERALL_FIELD_ll
!
IMPLICIT NONE
!
!
!* 0.1 declarations of arguments
!
-!
-REAL, INTENT(IN) :: PTSTEP ! time step
-REAL, INTENT(IN) :: PSEGLEN ! segment length
-INTEGER, INTENT(IN) :: KRR ! number of moist variables
-INTEGER, INTENT(IN) :: KSV ! number of scalar variables
-LOGICAL, INTENT(IN) :: OUSETKE ! flag to use tke
-REAL, INTENT(IN) :: PLATOR ! latitude of origine point
-REAL, INTENT(IN) :: PLONOR ! longitude of origine point
-INTEGER, INTENT(IN) :: KMI ! MODEL NUMBER
+! NONE
!
!-------------------------------------------------------------------------------
!
! 0.2 declaration of local variables
!
-INTEGER :: ISTORE ! number of storage instants
-INTEGER :: IIU_ll,IJU_ll,IRESP
+INTEGER :: IERR
+INTEGER :: IIU
+INTEGER :: IJU
+INTEGER :: ISTORE ! number of storage instants
INTEGER :: JI
+LOGICAL :: GALLOCX, GALLOCY
+LOGICAL :: GINSIDE ! True if station is inside physical domain of model
+LOGICAL :: GPRESENT ! True if station is present on the current process
+REAL :: ZXHATM_PHYS_MIN, ZYHATM_PHYS_MIN ! Minimum X coordinate of mass points in the physical domain
+REAL :: ZXHATM_PHYS_MAX, ZYHATM_PHYS_MAX ! Minimum X coordinate of mass points in the physical domain
+REAL, DIMENSION(SIZE(XXHAT)) :: ZXHATM ! mass point coordinates
+REAL, DIMENSION(SIZE(XYHAT)) :: ZYHATM ! mass point coordinates
+REAL, DIMENSION(:), POINTER :: ZXHAT_GLOB
+REAL, DIMENSION(:), POINTER :: ZYHAT_GLOB
+TYPE(TSTATIONDATA) :: TZSTATION
!
!----------------------------------------------------------------------------
-!
-!* 1. Default values
-! --------------
-!
-NUMBSTAT = 0
-TSTATIONS_TIME%XTSTEP = XTSTEP
-!
-!
-!* 3. Stations initialization
-! -----------------------
-!
-IF (CFILE_STAT=="NO_INPUT_CSV") THEN
-!
-!* 1. Namelist
-! --------
- NUMBSTAT = NNUMB_STAT
- IF (NUMBSTAT > 0) THEN
- ALLOCATE( TSTATIONS(NUMBSTAT) )
+TSTATIONS_TIME%XTSTEP = XSTEP_STAT
- IF (LCARTESIAN) THEN
- DO JI=1,NUMBSTAT
- TSTATIONS(JI)%XX = XX_STAT(JI)
- TSTATIONS(JI)%XY = XY_STAT(JI)
- TSTATIONS(JI)%XZ = XZ_STAT(JI)
- TSTATIONS(JI)%CNAME = CNAME_STAT(JI)
- TSTATIONS(JI)%CTYPE = CTYPE_STAT(JI)
- END DO
- ELSE
- DO JI=1,NUMBSTAT
- TSTATIONS(JI)%XLAT = XLAT_STAT(JI)
- TSTATIONS(JI)%XLON = XLON_STAT(JI)
- TSTATIONS(JI)%XZ = XZ_STAT(JI)
- TSTATIONS(JI)%CNAME = CNAME_STAT(JI)
- TSTATIONS(JI)%CTYPE = CTYPE_STAT(JI)
- END DO
- END IF
- END IF
-ELSE
-!
-!* 2. CSV DATA
-!
- CALL READ_CSV_STATION( CFILE_STAT, TSTATIONS, LCARTESIAN )
-END IF
+if ( tstations_time%xtstep < xtstep ) then
+ call Print_msg( NVERB_WARNING, 'GEN', 'INI_SURFSTATION_n', 'Timestep for stations was smaller than model timestep' )
+ tstations_time%xtstep = xtstep
+end if
-TSTATIONS_TIME%XTSTEP = XSTEP_STAT
+ISTORE = NINT ( ( XSEGLEN - DYN_MODEL(1)%XTSTEP ) / TSTATIONS_TIME%XTSTEP ) + 1
-LSTATION = (NUMBSTAT>0)
-!
-!----------------------------------------------------------------------------
+allocate( tstations_time%tpdates(istore) )
!
-!* 4. Allocations of storage arrays
-! -----------------------------
+! Prepare positioning data
!
-IF(NUMBSTAT>0) THEN
- CALL ALLOCATE_STATION_n()
- IF (.NOT. LCARTESIAN) CALL INI_INTERP_STATION_n()
-ENDIF
-!----------------------------------------------------------------------------
+IF ( CFILE_STAT /= "NO_INPUT_CSV" .OR. NNUMB_STAT > 0 ) THEN
+ IIU = SIZE( XXHAT )
+ IJU = SIZE( XYHAT )
+
+ ! Get global XHAT and YHAT (needed by STATION_POSITION)
+ CALL ALLOCBUFFER_ll( ZXHAT_GLOB, XXHAT, 'XX', GALLOCX )
+ CALL ALLOCBUFFER_ll( ZYHAT_GLOB, XYHAT, 'YY', GALLOCY )
+ CALL GATHERALL_FIELD_ll( 'XX', XXHAT, ZXHAT_GLOB, IERR )
+ CALL GATHERALL_FIELD_ll( 'YY', XYHAT, ZYHAT_GLOB, IERR )
+
+ ! 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)
+
+ ZXHATM_PHYS_MIN = 0.5 * ( ZXHAT_GLOB(1+JPHEXT) + ZXHAT_GLOB(2+JPHEXT) )
+ ZXHATM_PHYS_MAX = 0.5 * ( ZXHAT_GLOB(UBOUND(ZXHAT_GLOB,1)-JPHEXT) + ZXHAT_GLOB(UBOUND(ZXHAT_GLOB,1)-JPHEXT+1) )
+ ZYHATM_PHYS_MIN = 0.5 * ( ZYHAT_GLOB(1+JPHEXT) + ZYHAT_GLOB(2+JPHEXT) )
+ ZYHATM_PHYS_MAX = 0.5 * ( ZYHAT_GLOB(UBOUND(ZYHAT_GLOB,1)-JPHEXT) + ZYHAT_GLOB(UBOUND(ZYHAT_GLOB,1)-JPHEXT+1) )
+END IF
!
-CONTAINS
+! Stations initialization
!
-!----------------------------------------------------------------------------
-SUBROUTINE ALLOCATE_STATION_n()
+NUMBSTAT_LOC = 0
-INTEGER :: JI
-
-if ( tstations_time%xtstep < xtstep ) then
- call Print_msg( NVERB_WARNING, 'GEN', 'INI_SURFSTATION_n', 'Timestep for stations was smaller than model timestep' )
- tstations_time%xtstep = xtstep
-end if
+IF (CFILE_STAT=="NO_INPUT_CSV") THEN
+ ! Treat namelist
+ NUMBSTAT = 0
+ IF ( NNUMB_STAT > 0 ) THEN
+ DO JI = 1, NNUMB_STAT
+ IF ( LCARTESIAN ) THEN
+ TZSTATION%XX = XX_STAT(JI)
+ TZSTATION%XY = XY_STAT(JI)
+ ELSE
+ TZSTATION%XLAT = XLAT_STAT(JI)
+ TZSTATION%XLON = XLON_STAT(JI)
+ CALL STATION_INI_INTERP( TZSTATION )
+ END IF
+ TZSTATION%XZ = XZ_STAT(JI)
+ TZSTATION%CNAME = CNAME_STAT(JI)
-ISTORE = NINT ( ( PSEGLEN - DYN_MODEL(1)%XTSTEP ) / TSTATIONS_TIME%XTSTEP ) + 1
+ CALL STATION_POSITION( TZSTATION, ZXHAT_GLOB, ZYHAT_GLOB, ZXHATM, ZYHATM, &
+ ZXHATM_PHYS_MIN, ZXHATM_PHYS_MAX, ZYHATM_PHYS_MIN, ZYHATM_PHYS_MAX, &
+ GINSIDE, GPRESENT )
-allocate( tstations_time%tpdates(istore) )
+ IF ( GINSIDE ) THEN
+ NUMBSTAT = NUMBSTAT + 1
+ TZSTATION%NID = NUMBSTAT
+ END IF
-DO JI = 1, NUMBSTAT
- ALLOCATE(TSTATIONS(JI)%XZON (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XMER (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XW (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XP (ISTORE))
- IF (OUSETKE) THEN
- ALLOCATE(TSTATIONS(JI)%XTKE (ISTORE))
- ELSE
- ALLOCATE(TSTATIONS(JI)%XTKE (0))
+ IF ( GPRESENT ) CALL STATION_ADD( TZSTATION )
+ END DO
END IF
- ALLOCATE(TSTATIONS(JI)%XTH (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XR (ISTORE,KRR))
- ALLOCATE(TSTATIONS(JI)%XSV (ISTORE,KSV))
- ALLOCATE(TSTATIONS(JI)%XTSRAD (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XT2M (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XQ2M (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XHU2M (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XZON10M(ISTORE))
- ALLOCATE(TSTATIONS(JI)%XMER10M(ISTORE))
- ALLOCATE(TSTATIONS(JI)%XRN (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XH (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XLE (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XLEI (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XGFLUX (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XSWD (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XSWU (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XLWD (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XLWU (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XSWDIR (ISTORE))
- ALLOCATE(TSTATIONS(JI)%XSWDIFF(ISTORE))
- ALLOCATE(TSTATIONS(JI)%XDSTAOD(ISTORE))
- ALLOCATE(TSTATIONS(JI)%XSFCO2 (ISTORE))
+ELSE
+ !Treat CSV datafile
+ CALL READ_CSV_STATION( CFILE_STAT, ZXHAT_GLOB, ZYHAT_GLOB, ZXHATM, ZYHATM, &
+ ZXHATM_PHYS_MIN, ZXHATM_PHYS_MAX,ZYHATM_PHYS_MIN, ZYHATM_PHYS_MAX )
+END IF
- TSTATIONS(JI)%XZON(:) = XUNDEF
- TSTATIONS(JI)%XMER(:) = XUNDEF
- TSTATIONS(JI)%XW(:) = XUNDEF
- TSTATIONS(JI)%XP(:) = XUNDEF
- TSTATIONS(JI)%XTKE(:) = XUNDEF
- TSTATIONS(JI)%XTH(:) = XUNDEF
- TSTATIONS(JI)%XR(:,:) = XUNDEF
- TSTATIONS(JI)%XSV(:,:) = XUNDEF
- TSTATIONS(JI)%XTSRAD(:) = XUNDEF
- TSTATIONS(JI)%XT2M(:) = XUNDEF
- TSTATIONS(JI)%XQ2M(:) = XUNDEF
- TSTATIONS(JI)%XHU2M(:) = XUNDEF
- TSTATIONS(JI)%XZON10M(:) = XUNDEF
- TSTATIONS(JI)%XMER10M(:) = XUNDEF
- TSTATIONS(JI)%XRN(:) = XUNDEF
- TSTATIONS(JI)%XH(:) = XUNDEF
- TSTATIONS(JI)%XLE(:) = XUNDEF
- TSTATIONS(JI)%XLEI(:) = XUNDEF
- TSTATIONS(JI)%XGFLUX(:) = XUNDEF
- TSTATIONS(JI)%XSWD(:) = XUNDEF
- TSTATIONS(JI)%XSWU(:) = XUNDEF
- TSTATIONS(JI)%XLWD(:) = XUNDEF
- TSTATIONS(JI)%XLWU(:) = XUNDEF
- TSTATIONS(JI)%XSWDIR(:) = XUNDEF
- TSTATIONS(JI)%XSWDIFF(:) = XUNDEF
- TSTATIONS(JI)%XDSTAOD(:) = XUNDEF
- TSTATIONS(JI)%XSFCO2(:) = XUNDEF
-END DO
+LSTATION = ( NUMBSTAT > 0 )
-END SUBROUTINE ALLOCATE_STATION_n
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-SUBROUTINE INI_INTERP_STATION_n()
+DO JI = 1, NUMBSTAT_LOC
+ CALL STATION_ALLOCATE( TSTATIONS(JI), ISTORE )
+END DO
!
-USE MODE_STATION_TOOLS, ONLY: STATION_POSITION
-
-INTEGER :: JII
-INTEGER :: IIU, IJU
+! Clean positioning data
!
-IF ( ALL(TSTATIONS(:)%XLAT /= XUNDEF) .AND. ALL(TSTATIONS(:)%XLON /= XUNDEF) ) THEN
- DO JII = 1, NUMBSTAT
- CALL GET_DIM_EXT_ll ('B',IIU,IJU)
- 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'
- CMNHMSG(1) = 'either LATitude or LONgitude segment'
- CMNHMSG(1) = 'or I and J segment'
- CMNHMSG(1) = 'definition is not complete.'
- CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'INI_SURFSTATION_n' )
+IF ( CFILE_STAT /= "NO_INPUT_CSV" .OR. NNUMB_STAT > 0 ) THEN
+ IF ( GALLOCX ) DEALLOCATE( ZXHAT_GLOB )
+ IF ( GALLOCY ) DEALLOCATE( ZYHAT_GLOB )
END IF
-END SUBROUTINE INI_INTERP_STATION_n
!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-!
+
END SUBROUTINE INI_SURFSTATION_n
diff --git a/src/MNH/modd_allstationn.f90 b/src/MNH/modd_allstationn.f90
index 933c16571..4229177e4 100644
--- a/src/MNH/modd_allstationn.f90
+++ b/src/MNH/modd_allstationn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2021-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2021-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.
@@ -29,6 +29,7 @@
!! MODIFICATIONS
!! -------------
!! Original 01/06/21
+! P. Wautelet 04/2022: restructure stations for better performance, reduce memory usage and correct some problems/bugs
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -36,11 +37,16 @@
!
!
USE MODD_PARAMETERS, ONLY: JPMODELMAX
-USE MODD_STATION_n
-USE MODD_TYPE_STATION
IMPLICIT NONE
+PRIVATE
+
+PUBLIC :: NNUMB_STAT, XSTEP_STAT, XX_STAT, XY_STAT, XLAT_STAT, XLON_STAT, XZ_STAT
+PUBLIC :: CNAME_STAT, CFILE_STAT, LDIAG_SURFRAD
+
+PUBLIC :: ALLSTATION_GOTO_MODEL
+
TYPE ALLSTATION_t
!
!-------------------------------------------------------------------------------------------
@@ -48,7 +54,7 @@ TYPE ALLSTATION_t
!
INTEGER :: NNUMB_STAT !Number of stations as defined in namelist
REAL, DIMENSION(100) :: XX_STAT, XY_STAT, XZ_STAT, XLAT_STAT, XLON_STAT
- CHARACTER(LEN=7), DIMENSION(100) :: CNAME_STAT, CTYPE_STAT
+ CHARACTER(LEN=7), DIMENSION(100) :: CNAME_STAT
CHARACTER(LEN=20) :: CFILE_STAT
REAL :: XSTEP_STAT
LOGICAL :: LDIAG_SURFRAD
@@ -66,7 +72,6 @@ REAL, DIMENSION(:), POINTER :: XLAT_STAT=>NULL()
REAL, DIMENSION(:), POINTER :: XLON_STAT=>NULL()
REAL, DIMENSION(:), POINTER :: XZ_STAT=>NULL()
CHARACTER (LEN=7),DIMENSION(:), POINTER :: CNAME_STAT=>NULL()
-CHARACTER (LEN=7),DIMENSION(:), POINTER :: CTYPE_STAT=>NULL()
CHARACTER (LEN=20),POINTER :: CFILE_STAT=>NULL()
LOGICAL, POINTER :: LDIAG_SURFRAD=>NULL()
CONTAINS
@@ -86,7 +91,6 @@ XZ_STAT =>ALLSTATION_MODEL(KTO)%XZ_STAT
XLAT_STAT =>ALLSTATION_MODEL(KTO)%XLAT_STAT
XLON_STAT =>ALLSTATION_MODEL(KTO)%XLON_STAT
CNAME_STAT =>ALLSTATION_MODEL(KTO)%CNAME_STAT
-CTYPE_STAT =>ALLSTATION_MODEL(KTO)%CTYPE_STAT
CFILE_STAT =>ALLSTATION_MODEL(KTO)%CFILE_STAT
LDIAG_SURFRAD =>ALLSTATION_MODEL(KTO)%LDIAG_SURFRAD
END SUBROUTINE ALLSTATION_GOTO_MODEL
diff --git a/src/MNH/modd_stationn.f90 b/src/MNH/modd_stationn.f90
index ff5351dee..6017c0f98 100644
--- a/src/MNH/modd_stationn.f90
+++ b/src/MNH/modd_stationn.f90
@@ -29,25 +29,31 @@
!! MODIFICATIONS
!! -------------
!! Original 15/01/02
-! P. Wautelet 07/04/2022: rewrite types for stations
+! P. Wautelet 04/2022: restructure stations for better performance, reduce memory usage and correct some problems/bugs
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
!
-USE MODD_PARAMETERS, ONLY: JPMODELMAX
-USE MODD_TYPE_STATION
+USE MODD_PARAMETERS, ONLY: JPMODELMAX
+USE MODD_TYPE_STATION, ONLY: TSTATIONDATA, TSTATIONTIME
IMPLICIT NONE
+PRIVATE
+
+PUBLIC :: LSTATION, NUMBSTAT, NUMBSTAT_LOC, TSTATIONS_TIME, TSTATIONS
+
+PUBLIC :: STATION_GOTO_MODEL
+
TYPE STATION_t
!
!-------------------------------------------------------------------------------------------
!
LOGICAL :: LSTATION ! flag to use stations
INTEGER :: NUMBSTAT ! number of stations
- LOGICAL :: LSTATLAT ! positioning in lat/lon
+ INTEGER :: NUMBSTAT_LOC = 0 ! number of stations on this process
!
TYPE(TSTATIONTIME) :: TSTATIONS_TIME
TYPE(TSTATIONDATA), DIMENSION(:), POINTER :: TSTATIONS ! characteristics and records of the stations
@@ -58,7 +64,7 @@ TYPE(STATION_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: STATION_MODEL
LOGICAL, POINTER :: LSTATION=>NULL()
INTEGER, POINTER :: NUMBSTAT=>NULL()
-LOGICAL, POINTER :: LSTATLAT=>NULL()
+INTEGER, POINTER :: NUMBSTAT_LOC=>NULL()
TYPE(TSTATIONTIME), POINTER :: TSTATIONS_TIME => NULL()
TYPE(TSTATIONDATA), DIMENSION(:), POINTER :: TSTATIONS => NULL()
@@ -73,7 +79,7 @@ STATION_MODEL(KFROM)%TSTATIONS => TSTATIONS
! Current model is set to model KTO
LSTATION => STATION_MODEL(KTO)%LSTATION
NUMBSTAT => STATION_MODEL(KTO)%NUMBSTAT
-LSTATLAT => STATION_MODEL(KTO)%LSTATLAT
+NUMBSTAT_LOC => STATION_MODEL(KTO)%NUMBSTAT_LOC
TSTATIONS_TIME => STATION_MODEL(KTO)%TSTATIONS_TIME
TSTATIONS => STATION_MODEL(KTO)%TSTATIONS
diff --git a/src/MNH/modd_type_station.f90 b/src/MNH/modd_type_station.f90
index ac0134eab..8d34883fa 100644
--- a/src/MNH/modd_type_station.f90
+++ b/src/MNH/modd_type_station.f90
@@ -30,14 +30,14 @@
!! -------------
!! Original 15/01/02
! P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
-! P. Wautelet 07/04/2022: rewrite types for stations
+! P. Wautelet 04/2022: restructure stations for better performance, reduce memory usage and correct some problems/bugs
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
use modd_type_date, only: date_time
-use modd_parameters, only: NNEGUNDEF, NUNDEF, XUNDEF
+use modd_parameters, only: NNEGUNDEF, XUNDEF
implicit none
@@ -54,10 +54,10 @@ END TYPE TSTATIONTIME
TYPE TSTATIONDATA
! Type to store all the data of 1 station
+
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
+
+INTEGER :: NID = 0 ! Global identification number of the station (from 1 to total number of stations of the model)
REAL :: XX = XUNDEF ! X(n)
REAL :: XY = XUNDEF ! Y(n)
@@ -78,7 +78,7 @@ 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
+INTEGER :: NK = NNEGUNDEF ! Model level for altitude comparisons
REAL, DIMENSION(:), ALLOCATABLE :: XZON ! zonal wind(n)
REAL, DIMENSION(:), ALLOCATABLE :: XMER ! meridian wind(n)
diff --git a/src/MNH/modn_stationn.f90 b/src/MNH/modn_stationn.f90
index f388061e7..ab6013ee2 100644
--- a/src/MNH/modn_stationn.f90
+++ b/src/MNH/modn_stationn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2020-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2020-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.
@@ -21,6 +21,7 @@
!! MODIFICATIONS
!! -------------
!! Original 10/03/20
+! P. Wautelet 04/2022: restructure stations for better performance, reduce memory usage and correct some problems/bugs
!!
!! IMPLICIT ARGUMENTS
!! ------------------
@@ -34,7 +35,6 @@ USE MODD_ALLSTATION_n, ONLY:&
XLON_STAT_n =>XLON_STAT ,&
XZ_STAT_n =>XZ_STAT ,&
CNAME_STAT_n =>CNAME_STAT ,&
- CTYPE_STAT_n =>CTYPE_STAT ,&
CFILE_STAT_n =>CFILE_STAT ,&
LDIAG_SURFRAD_n =>LDIAG_SURFRAD
!!
@@ -46,7 +46,7 @@ IMPLICIT NONE
INTEGER ,SAVE:: NNUMB_STAT
REAL ,SAVE:: XSTEP_STAT
REAL, DIMENSION(100) ,SAVE:: XX_STAT, XY_STAT, XZ_STAT, XLAT_STAT, XLON_STAT
-CHARACTER (LEN=7), DIMENSION(100),SAVE:: CNAME_STAT, CTYPE_STAT
+CHARACTER (LEN=7), DIMENSION(100),SAVE:: CNAME_STAT
CHARACTER (LEN=20) ,SAVE:: CFILE_STAT !filename
LOGICAL ,SAVE:: LDIAG_SURFRAD
@@ -54,7 +54,7 @@ NAMELIST /NAM_STATIONn/ &
NNUMB_STAT, XSTEP_STAT, &
XX_STAT,XY_STAT,XZ_STAT,&
XLON_STAT,XLAT_STAT,&
- CNAME_STAT,CTYPE_STAT,&
+ CNAME_STAT,&
CFILE_STAT,LDIAG_SURFRAD
!
@@ -69,7 +69,6 @@ SUBROUTINE INIT_NAM_STATIONn
XLON_STAT = XLON_STAT_n
XZ_STAT = XZ_STAT_n
CNAME_STAT = CNAME_STAT_n
- CTYPE_STAT = CTYPE_STAT_n
CFILE_STAT = CFILE_STAT_n
LDIAG_SURFRAD= LDIAG_SURFRAD_n
END SUBROUTINE INIT_NAM_STATIONn
@@ -83,7 +82,6 @@ SUBROUTINE UPDATE_NAM_STATIONn
XLON_STAT_n = XLON_STAT
XZ_STAT_n = XZ_STAT
CNAME_STAT_n = CNAME_STAT
- CTYPE_STAT_n = CTYPE_STAT
CFILE_STAT_n = CFILE_STAT
LDIAG_SURFRAD_n= LDIAG_SURFRAD
END SUBROUTINE UPDATE_NAM_STATIONn
diff --git a/src/MNH/station_reader.f90 b/src/MNH/station_reader.f90
index 3e7741e8b..1a72bfe1a 100644
--- a/src/MNH/station_reader.f90
+++ b/src/MNH/station_reader.f90
@@ -4,21 +4,18 @@
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
! #######################
- MODULE MODI_STATION_READER
+ MODULE MODE_STATION_READER
! #######################
-!
-INTERFACE
-!
-SUBROUTINE READ_CSV_STATION( HFILE, TPSTATIONS, OCARTESIAN )
-USE MODD_STATION_n
-CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
-TYPE(TSTATIONDATA), DIMENSION(:), POINTER, INTENT(INOUT) :: TPSTATIONS
-LOGICAL, INTENT(IN) :: OCARTESIAN
-END SUBROUTINE READ_CSV_STATION
-!
-END INTERFACE
-!
-END MODULE MODI_STATION_READER
+
+IMPLICIT NONE
+
+PRIVATE
+
+PUBLIC :: READ_CSV_STATION
+
+INTEGER, PARAMETER :: NMAXLINELGT = 400
+
+CONTAINS
!-------------------------------------------------------------------
!
!!**** *READ_CSV_STATION* -
@@ -34,85 +31,85 @@ END MODULE MODI_STATION_READER
!! MODIFICATIONS
!! -------------
!! 03/2020 Original
-! P. Wautelet 07/04/2022: rewrite types for stations
+! P. Wautelet 04/2022: restructure stations for better performance, reduce memory usage and correct some problems/bugs
!---------------------------------------------------------------
!
-!#########################################################
-SUBROUTINE READ_CSV_STATION( HFILE, TPSTATIONS, OCARTESIAN )
-USE MODD_ALLSTATION_n
-USE MODD_STATION_n
-USE MODD_PARAMETERS
-USE MODD_TYPE_STATION
-USE MODI_INI_SURFSTATION_n
+!###############################################################################################
+SUBROUTINE READ_CSV_STATION( HFILE, PXHAT_GLOB, PYHAT_GLOB, PXHATM, PYHATM, &
+ PXHATM_PHYS_MIN, PXHATM_PHYS_MAX,PYHATM_PHYS_MIN, PYHATM_PHYS_MAX )
+!###############################################################################################
+USE MODD_CONF, ONLY: LCARTESIAN
+USE MODD_STATION_n, ONLY: NUMBSTAT
+USE MODD_TYPE_STATION, ONLY: TSTATIONDATA
+
+USE MODE_MSG
+USE MODE_STATION_TOOLS, ONLY: STATION_ADD, STATION_INI_INTERP, STATION_POSITION
+
+CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+REAL, DIMENSION(:), INTENT(IN) :: PXHAT_GLOB
+REAL, DIMENSION(:), INTENT(IN) :: PYHAT_GLOB
+REAL, DIMENSION(:), INTENT(IN) :: PXHATM ! mass point coordinates
+REAL, DIMENSION(:), INTENT(IN) :: PYHATM ! mass point coordinates
+REAL, INTENT(IN) :: PXHATM_PHYS_MIN, PYHATM_PHYS_MIN ! Minimum X coordinate of mass points in the physical domain
+REAL, INTENT(IN) :: PXHATM_PHYS_MAX, PYHATM_PHYS_MAX ! Minimum X coordinate of mass points in the physical domain
!
-CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
-TYPE(TSTATIONDATA), DIMENSION(:), POINTER, INTENT(INOUT) :: TPSTATIONS
-LOGICAL, INTENT(IN) :: OCARTESIAN
-!
-INTEGER :: INBLINE ! Nb of line in csv file
-!
-CHARACTER(LEN=80) :: YERROR
-CHARACTER(LEN=400) :: YSTRING
-INTEGER :: ILU ! logical unit of the file
-!
+CHARACTER(LEN=NMAXLINELGT) :: YSTRING
+INTEGER :: ILU ! logical unit of the file
+INTEGER :: INBLINE ! Nb of lines in csv file
+INTEGER :: JI
+LOGICAL :: GINSIDE ! True if station is inside physical domain of model
+LOGICAL :: GPRESENT ! True if station is present on the current process
+TYPE(TSTATIONDATA) :: TZSTATION
+
+INBLINE = 0 !Number of stations found in the file
+NUMBSTAT = 0 !Number of stations found in the file AND inside the model domain
! Open file
-OPEN(NEWUNIT=ILU,FILE=HFILE, FORM='formatted')
-! Count lines
-REWIND(ILU)
-INBLINE=0
+OPEN( NEWUNIT = ILU, FILE = HFILE, FORM = 'formatted' )
+
+READ( ILU, END = 101, FMT = '(A)' ) YSTRING ! Reading of header (skip it)
+
DO
- READ(ILU,END=101,FMT='(A400)') YSTRING
-!* analyses if the record has been written in French convention
- CALL FRENCH_TO_ENGLISH(YSTRING) ! analyse de convention fr ou eng
- IF (LEN_TRIM(YSTRING) > 0) THEN
+ ! Read station coordinates
+ READ( ILU, END = 101, FMT = '(A)' ) YSTRING
+
+ ! Check if record is written in French convention
+ CALL FRENCH_TO_ENGLISH( YSTRING )
+
+ IF ( LCARTESIAN ) THEN
+ READ( YSTRING, * ) TZSTATION%CNAME, TZSTATION%XX, TZSTATION%XY, TZSTATION%XZ
+ ELSE
+ READ( YSTRING, * ) TZSTATION%CNAME, TZSTATION%XLAT, TZSTATION%XLON, TZSTATION%XZ
+ END IF
+
+ IF ( .NOT. LCARTESIAN ) CALL STATION_INI_INTERP( TZSTATION )
+ CALL STATION_POSITION( TZSTATION, PXHAT_GLOB, PYHAT_GLOB, PXHATM, PYHATM, &
+ PXHATM_PHYS_MIN, PXHATM_PHYS_MAX, PYHATM_PHYS_MIN, PYHATM_PHYS_MAX, &
+ GINSIDE, GPRESENT )
+
+ IF ( GINSIDE ) THEN
+ NUMBSTAT = NUMBSTAT + 1
+ TZSTATION%NID = NUMBSTAT
+ END IF
+
+ IF ( GPRESENT ) CALL STATION_ADD( TZSTATION )
+
INBLINE = INBLINE + 1
- END IF
END DO
-!
+
101 CONTINUE
- IF (INBLINE == 0) THEN
- YERROR = 'Data not found in file : '//TRIM(HFILE)
- PRINT*, YERROR
- ELSE
- ! Save number of stations
- NUMBSTAT = INBLINE - 1
-
- ALLOCATE( TPSTATIONS(NUMBSTAT) )
-
- ! New reading
- REWIND(ILU)
- READ(ILU,FMT='(A400)') YSTRING ! Reading of header
- !
- ! Save the data
- IF (OCARTESIAN) THEN
- INBLINE = 1
- DO INBLINE=1, NUMBSTAT
- READ(ILU,FMT='(A400)') YSTRING
- READ(YSTRING,*) TPSTATIONS(INBLINE)%CNAME, & !TPSTATIONS(INBLINE)%CTYPE,&
- TPSTATIONS(INBLINE)%XX, TPSTATIONS(INBLINE)%XY, TPSTATIONS(INBLINE)%XZ
- END DO
- REWIND(ILU)
- CLOSE(ILU)
- RETURN
- ELSE
- INBLINE = 1
- DO INBLINE=1, NUMBSTAT
- READ(ILU,FMT='(A400)') YSTRING
- READ(YSTRING,*) TPSTATIONS(INBLINE)%CNAME, & !TPSTATIONS(INBLINE)%CTYPE,&
- TPSTATIONS(INBLINE)%XLAT, TPSTATIONS(INBLINE)%XLON, TPSTATIONS(INBLINE)%XZ
- END DO
- REWIND(ILU)
- CLOSE(ILU)
- RETURN
- END IF
- END IF
-!
+
+CLOSE( ILU )
+
+IF ( INBLINE == 0 ) CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'READ_CSV_STATION', 'Data not found in file ' // TRIM( HFILE ) )
+
END SUBROUTINE READ_CSV_STATION
+
!#########################################################
SUBROUTINE FRENCH_TO_ENGLISH(HSTRING)
-CHARACTER(LEN=400), INTENT(INOUT) :: HSTRING ! csv record
+CHARACTER(LEN=NMAXLINELGT), INTENT(INOUT) :: HSTRING ! csv record
+
INTEGER :: JL
LOGICAL :: GFRENCH
!
@@ -120,13 +117,13 @@ GFRENCH = .FALSE.
!* analyses if the record has been written in French convention
! French convention (separator is ; decimal symbol is ,)
! or English convention (separator is , decimal symbol is .)
-DO JL=1,400
+DO JL = 1, NMAXLINELGT
IF (HSTRING(JL:JL)==';') GFRENCH=.TRUE.
END DO
!
! If French convention is used in the file, transforms it in English convention
IF (GFRENCH) THEN
- DO JL=1,400
+ DO JL = 1, NMAXLINELGT
IF (HSTRING(JL:JL)==',') HSTRING(JL:JL)='.'
IF (HSTRING(JL:JL)==';') HSTRING(JL:JL)=','
END DO
@@ -134,3 +131,4 @@ END IF
!
END SUBROUTINE FRENCH_TO_ENGLISH
+END MODULE MODE_STATION_READER
diff --git a/src/MNH/station_tools.f90 b/src/MNH/station_tools.f90
index 1aa02f68e..10b3c1ee3 100644
--- a/src/MNH/station_tools.f90
+++ b/src/MNH/station_tools.f90
@@ -1,9 +1,10 @@
-!MNH_LIC Copyright 2022-2022 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2002-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:
+! Authors:
+! Misc: some of the code was taken from older subroutines/functions for stations
! P. Wautelet 08/04/2022
!-----------------------------------------------------------------
! Modifications:
@@ -12,84 +13,381 @@
MODULE MODE_STATION_TOOLS
! ##################
+USE MODD_TYPE_STATION, ONLY: TSTATIONDATA
+
IMPLICIT NONE
PRIVATE
+PUBLIC :: STATION_ALLOCATE
+PUBLIC :: STATION_INI_INTERP
PUBLIC :: STATION_POSITION
+PUBLIC :: STATION_ADD
+PUBLIC :: STATION_INTERP_2D, STATION_INTERP_2D_U, STATION_INTERP_2D_V
CONTAINS
-! ######################################
-SUBROUTINE STATION_POSITION( TPSTATION )
-! ######################################
+! ##############################################
+SUBROUTINE STATION_ALLOCATE( TPSTATION, KSTORE )
+! ##############################################
+
+ USE MODD_ALLSTATION_n, ONLY: LDIAG_SURFRAD
+ USE MODD_CONF_n, ONLY: NRR
+ USE MODD_NSV, ONLY: NSV
+ USE MODD_PARAMETERS, ONLY: XUNDEF
+ USE MODD_PARAM_n, ONLY: CRAD, CTURB
+
+ IMPLICIT NONE
+
+ TYPE(TSTATIONDATA), INTENT(INOUT) :: TPSTATION
+ INTEGER, INTENT(IN) :: KSTORE ! number of moments to store
+
+ ALLOCATE( TPSTATION%XZON(KSTORE) )
+ ALLOCATE( TPSTATION%XMER(KSTORE) )
+ ALLOCATE( TPSTATION%XW (KSTORE) )
+ ALLOCATE( TPSTATION%XP (KSTORE) )
+ IF ( CTURB == 'TKEL' ) THEN
+ ALLOCATE( TPSTATION%XTKE(KSTORE) )
+ ELSE
+ ALLOCATE( TPSTATION%XTKE(0) )
+ END IF
+ ALLOCATE( TPSTATION%XTH(KSTORE) )
+ ALLOCATE( TPSTATION%XR (KSTORE, NRR) )
+ ALLOCATE( TPSTATION%XSV(KSTORE, NSV) )
+ IF ( CRAD /= 'NONE' ) THEN
+ ALLOCATE( TPSTATION%XTSRAD(KSTORE) )
+ ELSE
+ ALLOCATE( TPSTATION%XTSRAD(0) )
+ END IF
+ IF ( LDIAG_SURFRAD ) THEN
+ ALLOCATE( TPSTATION%XT2M (KSTORE) )
+ ALLOCATE( TPSTATION%XQ2M (KSTORE) )
+ ALLOCATE( TPSTATION%XHU2M (KSTORE) )
+ ALLOCATE( TPSTATION%XZON10M(KSTORE) )
+ ALLOCATE( TPSTATION%XMER10M(KSTORE) )
+ ALLOCATE( TPSTATION%XRN (KSTORE) )
+ ALLOCATE( TPSTATION%XH (KSTORE) )
+ ALLOCATE( TPSTATION%XLE (KSTORE) )
+ ALLOCATE( TPSTATION%XLEI (KSTORE) )
+ ALLOCATE( TPSTATION%XGFLUX (KSTORE) )
+ IF ( CRAD /= 'NONE' ) THEN
+ ALLOCATE( TPSTATION%XSWD (KSTORE) )
+ ALLOCATE( TPSTATION%XSWU (KSTORE) )
+ ALLOCATE( TPSTATION%XLWD (KSTORE) )
+ ALLOCATE( TPSTATION%XLWU (KSTORE) )
+ ALLOCATE( TPSTATION%XSWDIR (KSTORE) )
+ ALLOCATE( TPSTATION%XSWDIFF(KSTORE) )
+ ALLOCATE( TPSTATION%XDSTAOD(KSTORE) )
+ END IF
+ ALLOCATE( TPSTATION%XSFCO2(KSTORE) )
+ END IF
+
+ TPSTATION%XZON(:) = XUNDEF
+ TPSTATION%XMER(:) = XUNDEF
+ TPSTATION%XW(:) = XUNDEF
+ TPSTATION%XP(:) = XUNDEF
+ TPSTATION%XTKE(:) = XUNDEF
+ TPSTATION%XTH(:) = XUNDEF
+ TPSTATION%XR(:,:) = XUNDEF
+ TPSTATION%XSV(:,:) = XUNDEF
+ TPSTATION%XTSRAD(:) = XUNDEF
+ IF ( LDIAG_SURFRAD ) THEN
+ TPSTATION%XT2M(:) = XUNDEF
+ TPSTATION%XQ2M(:) = XUNDEF
+ TPSTATION%XHU2M(:) = XUNDEF
+ TPSTATION%XZON10M(:) = XUNDEF
+ TPSTATION%XMER10M(:) = XUNDEF
+ TPSTATION%XRN(:) = XUNDEF
+ TPSTATION%XH(:) = XUNDEF
+ TPSTATION%XLE(:) = XUNDEF
+ TPSTATION%XLEI(:) = XUNDEF
+ TPSTATION%XGFLUX(:) = XUNDEF
+ IF ( CRAD /= 'NONE' ) THEN
+ TPSTATION%XSWD(:) = XUNDEF
+ TPSTATION%XSWU(:) = XUNDEF
+ TPSTATION%XLWD(:) = XUNDEF
+ TPSTATION%XLWU(:) = XUNDEF
+ TPSTATION%XSWDIR(:) = XUNDEF
+ TPSTATION%XSWDIFF(:) = XUNDEF
+ TPSTATION%XDSTAOD(:) = XUNDEF
+ END IF
+ TPSTATION%XSFCO2(:) = XUNDEF
+ END IF
+
+END SUBROUTINE STATION_ALLOCATE
+
+! ########################################
+SUBROUTINE STATION_INI_INTERP( TPSTATION )
+! ########################################
+
+ USE MODD_GRID, ONLY: XLATORI, XLONORI
+ USE MODD_PARAMETERS, ONLY: XUNDEF
+
+ USE MODE_GRIDPROJ, ONLY: SM_XYHAT
+ USE MODE_MSG
+
+ IMPLICIT NONE
+
+ TYPE(TSTATIONDATA), INTENT(INOUT) :: TPSTATION
+
+ IF ( TPSTATION%XLAT == XUNDEF .OR. TPSTATION%XLON == XUNDEF ) THEN
+ CMNHMSG(1) = 'Error in station position'
+ CMNHMSG(2) = 'either LATitude or LONgitude segment'
+ CMNHMSG(3) = 'or I and J segment'
+ CMNHMSG(4) = 'definition is not complete.'
+ CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'STATION_INI_INTERP' )
+ END IF
+
+ CALL SM_XYHAT( XLATORI, XLONORI, &
+ TPSTATION%XLAT, TPSTATION%XLON, &
+ TPSTATION%XX, TPSTATION%XY )
+
+END SUBROUTINE STATION_INI_INTERP
+
+! ###############################################################################################
+SUBROUTINE STATION_POSITION( TPSTATION, PXHAT_GLOB, PYHAT_GLOB, PXHATM, PYHATM, &
+ PXHATM_PHYS_MIN, PXHATM_PHYS_MAX,PYHATM_PHYS_MIN, PYHATM_PHYS_MAX, &
+ OINSIDE, OPRESENT )
+! ###############################################################################################
! Subroutine to determine the position of a station on the model grid
-! and set the useful coefficient for data interpolation
+! and set the useful coefficients for data interpolation
- USE MODD_CONF, ONLY: L1D
- USE MODD_GRID_n, ONLY: XXHAT, XYHAT
- USE MODD_TYPE_STATION, ONLY: TSTATIONDATA
+ USE MODD_CONF, ONLY: L1D
+ USE MODD_GRID_n, ONLY: XXHAT, XYHAT, XZZ
+ USE MODD_PARAMETERS, ONLY: JPHEXT, JPVEXT
- USE MODE_TOOLS_ll, ONLY: GET_INDICE_ll, LWEST_ll, LEAST_ll, LNORTH_ll, LSOUTH_ll
+ USE MODE_MSG
+ USE MODE_NEST_LL, ONLY: GET_MODEL_NUMBER_ll
+ USE MODE_TOOLS_ll, ONLY: GET_INDICE_ll
IMPLICIT NONE
TYPE(TSTATIONDATA), INTENT(INOUT) :: TPSTATION
+ REAL, DIMENSION(:), INTENT(IN) :: PXHAT_GLOB
+ REAL, DIMENSION(:), INTENT(IN) :: PYHAT_GLOB
+ REAL, DIMENSION(:), INTENT(IN) :: PXHATM ! mass point coordinates
+ REAL, DIMENSION(:), INTENT(IN) :: PYHATM ! mass point coordinates
+ REAL, INTENT(IN) :: PXHATM_PHYS_MIN, PYHATM_PHYS_MIN ! Minimum X coordinate of mass points in the physical domain
+ REAL, INTENT(IN) :: PXHATM_PHYS_MAX, PYHATM_PHYS_MAX ! Minimum X coordinate of mass points in the physical domain
+ LOGICAL, INTENT(OUT) :: OINSIDE ! True if station is inside physical domain of model
+ LOGICAL, INTENT(OUT) :: OPRESENT ! True if station is present on the current process
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 )
+ INTEGER :: IMI
+ INTEGER :: JK
+ REAL :: ZLOW, ZHIGH
- CALL GET_INDICE_ll (IIB, IJB, IIE, IJE )
+ OPRESENT = .FALSE.
+ OINSIDE = .FALSE.
- ! 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)
+ CALL GET_INDICE_ll( IIB, IJB, IIE, IJE )
- 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.
+ IF ( TPSTATION%XX >= PXHAT_GLOB(JPHEXT+1) .AND. TPSTATION%XX <= PXHAT_GLOB(UBOUND(PXHAT_GLOB,1)-JPHEXT+1) &
+ .AND. TPSTATION%XY >= PYHAT_GLOB(JPHEXT+1) .AND. TPSTATION%XY <= PYHAT_GLOB(UBOUND(PYHAT_GLOB,1)-JPHEXT+1) ) THEN
+ OINSIDE = .TRUE.
+ ELSE
+ CALL GET_MODEL_NUMBER_ll(IMI)
+ WRITE( CMNHMSG(1), "( 'station ', A, ' is outside of physical domain of model', I3 )" ) TRIM(TPSTATION%CNAME), IMI
+ CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'STATION_POSITION' )
+ END IF
! 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.
+ TPSTATION%NI_M = COUNT( PXHATM(:) <= TPSTATION%XX )
! Y position
TPSTATION%NJ_V = COUNT( XYHAT (:) <= TPSTATION%XY )
- TPSTATION%NJ_M = COUNT( ZYHATM(:) <= TPSTATION%XY )
+ TPSTATION%NJ_M = COUNT( PYHATM(:) <= 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
+ ! Position of station according to process
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.
+ .AND. TPSTATION%NJ_V >= IJB .AND. TPSTATION%NJ_V <= IJE ) OPRESENT = .TRUE.
+ IF ( L1D ) OPRESENT = .TRUE.
+
+ ! Check if station is too near of physical domain border (outside of physical domain for mass points)
+ IF ( OINSIDE .AND. .NOT. L1D ) THEN
+ IF ( TPSTATION%XX < PXHATM_PHYS_MIN .OR. TPSTATION%XX > PXHATM_PHYS_MAX &
+ .OR. TPSTATION%XY < PYHATM_PHYS_MIN .OR. TPSTATION%XY > PYHATM_PHYS_MAX ) THEN
+ CALL GET_MODEL_NUMBER_ll(IMI)
+ WRITE( CMNHMSG(1), "( 'station ', A, ' is outside of mass-points physical domain of model', I3 )" ) &
+ TRIM(TPSTATION%CNAME), IMI
+ CMNHMSG(2) = 'but is inside of flux-points physical domain.'
+ CMNHMSG(3) = 'Meaning: station is too close to the boundaries of physical domain.'
+ CMNHMSG(4) = '=> station disabled (not computed)'
+ CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'STATION_POSITION' )
+ OPRESENT = .FALSE.
+ OINSIDE = .FALSE.
+ END IF
+ END IF
! Computations only on correct process
- IF ( TPSTATION%LPRESENT .AND. .NOT. L1D ) THEN
+ IF ( OPRESENT .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) )
+ TPSTATION%XXMCOEF = ( TPSTATION%XX - PXHATM(TPSTATION%NI_M) ) / ( PXHATM(TPSTATION%NI_M+1) - PXHATM(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) )
+ TPSTATION%XYMCOEF = ( TPSTATION%XY - PYHATM(TPSTATION%NJ_M) ) / ( PYHATM(TPSTATION%NJ_M+1) - PYHATM(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
+ IF ( OPRESENT ) THEN
+ ! The closest K-level to the station altitude is chosen
+ JK = JPVEXT + 1
+ DO WHILE ( ( STATION_INTERP_2D( TPSTATION, XZZ(:,:,JK) ) - STATION_INTERP_2D( TPSTATION, XZZ(:,:,JPVEXT+1) ) ) < TPSTATION%XZ)
+ JK = JK + 1
+ END DO
+ ZLOW = STATION_INTERP_2D( TPSTATION, XZZ(:,:,JK-1) ) - STATION_INTERP_2D( TPSTATION, XZZ(:,:,JPVEXT+1) )
+ ZHIGH = STATION_INTERP_2D( TPSTATION, XZZ(:,:,JK ) ) - STATION_INTERP_2D( TPSTATION, XZZ(:,:,JPVEXT+1) )
+ !If the station is nearer from the lower level, select it
+ IF ( ( ZHIGH - TPSTATION%XZ ) > ( TPSTATION%XZ - ZLOW ) ) JK = JK - 1
+ TPSTATION%NK = JK
+ END IF
+
END SUBROUTINE STATION_POSITION
+! #################################
+SUBROUTINE STATION_ADD( TPSTATION )
+! #################################
+! Subroutine to add a station to the local list of stations
+ USE MODD_STATION_n, ONLY: NUMBSTAT_LOC, TSTATIONS
+
+ IMPLICIT NONE
+
+ TYPE(TSTATIONDATA), INTENT(IN) :: TPSTATION
+
+ INTEGER :: JS
+ TYPE(TSTATIONDATA), DIMENSION(:), POINTER :: TZSTATIONS
+
+ NUMBSTAT_LOC = NUMBSTAT_LOC + 1
+
+ ALLOCATE( TZSTATIONS( NUMBSTAT_LOC ) )
+ DO JS = 1, NUMBSTAT_LOC - 1
+ TZSTATIONS(JS) = TSTATIONS(JS)
+ END DO
+ TZSTATIONS(NUMBSTAT_LOC) = TPSTATION
+ IF ( ASSOCIATED( TSTATIONS ) ) DEALLOCATE( TSTATIONS ) !Can be done without memory leak because allocatable arrays were
+ !not yet allocated (will be done in STATION_ALLOCATE)
+ TSTATIONS => TZSTATIONS
+
+END SUBROUTINE STATION_ADD
+!
+! ######################################################
+FUNCTION STATION_INTERP_2D( TPSTATION, PA ) RESULT( PB )
+! ######################################################
+ USE MODD_CONF, ONLY: L1D
+
+ USE MODE_MSG
+
+ IMPLICIT NONE
+
+ TYPE(TSTATIONDATA), INTENT(IN) :: TPSTATION
+ REAL, DIMENSION(:,:), INTENT(IN) :: PA
+ REAL :: PB
+
+ INTEGER :: JI, JJ
+
+ IF ( SIZE( PA, 1 ) == 2 ) THEN
+ JI = 1
+ JJ = 1
+ ELSE IF ( L1D ) THEN
+ JI = 2
+ JJ = 2
+ ELSE
+ JI = TPSTATION%NI_M
+ JJ = TPSTATION%NJ_M
+ END IF
+
+ IF ( JI >= 1 .AND. JI < SIZE( PA, 1 ) &
+ .AND. JJ >= 1 .AND. JJ < SIZE( PA, 2 ) ) THEN
+ PB = (1.-TPSTATION%XXMCOEF) * (1.-TPSTATION%XYMCOEF) * PA(JI,JJ) + &
+ ( TPSTATION%XXMCOEF) * (1.-TPSTATION%XYMCOEF) * PA(JI+1,JJ) + &
+ (1.-TPSTATION%XXMCOEF) * ( TPSTATION%XYMCOEF) * PA(JI,JJ+1) + &
+ ( TPSTATION%XXMCOEF) * ( TPSTATION%XYMCOEF) * PA(JI+1,JJ+1)
+ ELSE
+ CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATION_INTERP_2D', 'value can not be interpolated' )
+ END IF
+
+END FUNCTION STATION_INTERP_2D
+
+! ########################################################
+FUNCTION STATION_INTERP_2D_U( TPSTATION, PA ) RESULT( PB )
+! ########################################################
+ USE MODD_CONF, ONLY: L1D
+
+ USE MODE_MSG
+
+ IMPLICIT NONE
+
+ TYPE(TSTATIONDATA), INTENT(IN) :: TPSTATION
+ REAL, DIMENSION(:,:), INTENT(IN) :: PA
+ REAL :: PB
+
+ INTEGER :: JI, JJ
+
+ IF ( SIZE( PA, 1 ) == 2 ) THEN
+ JI = 1
+ JJ = 1
+ ELSE IF ( L1D ) THEN
+ JI = 2
+ JJ = 2
+ ELSE
+ JI = TPSTATION%NI_U
+ JJ = TPSTATION%NJ_M
+ END IF
+
+ IF ( JI >= 1 .AND. JI < SIZE( PA, 1 ) &
+ .AND. JJ >= 1 .AND. JJ < SIZE( PA, 2 ) ) THEN
+ PB = (1.-TPSTATION%XXUCOEF) * (1.-TPSTATION%XYMCOEF) * PA(JI,JJ) + &
+ ( TPSTATION%XXUCOEF) * (1.-TPSTATION%XYMCOEF) * PA(JI+1,JJ) + &
+ (1.-TPSTATION%XXUCOEF) * ( TPSTATION%XYMCOEF) * PA(JI,JJ+1) + &
+ ( TPSTATION%XXUCOEF) * ( TPSTATION%XYMCOEF) * PA(JI+1,JJ+1)
+ ELSE
+ CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATION_INTERP_2D_U', 'value can not be interpolated' )
+ END IF
+
+END FUNCTION STATION_INTERP_2D_U
+
+! ########################################################
+FUNCTION STATION_INTERP_2D_V( TPSTATION, PA ) RESULT( PB )
+! ########################################################
+ USE MODD_CONF, ONLY: L1D
+
+ USE MODE_MSG
+
+ IMPLICIT NONE
+
+ TYPE(TSTATIONDATA), INTENT(IN) :: TPSTATION
+ REAL, DIMENSION(:,:), INTENT(IN) :: PA
+ REAL :: PB
+
+ INTEGER :: JI, JJ
+
+ IF ( SIZE( PA, 1 ) == 2 ) THEN
+ JI = 1
+ JJ = 1
+ ELSE IF ( L1D ) THEN
+ JI = 2
+ JJ = 2
+ ELSE
+ JI = TPSTATION%NI_M
+ JJ = TPSTATION%NJ_V
+ END IF
+
+ IF ( JI >= 1 .AND. JI < SIZE( PA, 1 ) &
+ .AND. JJ >= 1 .AND. JJ < SIZE( PA, 2 ) ) THEN
+ PB = (1.-TPSTATION%XXMCOEF) * (1.-TPSTATION%XYVCOEF) * PA(JI,JJ) + &
+ ( TPSTATION%XXMCOEF) * (1.-TPSTATION%XYVCOEF) * PA(JI+1,JJ) + &
+ (1.-TPSTATION%XXMCOEF) * ( TPSTATION%XYVCOEF) * PA(JI,JJ+1) + &
+ ( TPSTATION%XXMCOEF) * ( TPSTATION%XYVCOEF) * PA(JI+1,JJ+1)
+ ELSE
+ CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATION_INTERP_2D_V', 'value can not be interpolated' )
+ END IF
+
+END FUNCTION STATION_INTERP_2D_V
+
END MODULE MODE_STATION_TOOLS
diff --git a/src/MNH/stationn.f90 b/src/MNH/stationn.f90
index dcb38f6d3..0343d9da3 100644
--- a/src/MNH/stationn.f90
+++ b/src/MNH/stationn.f90
@@ -75,33 +75,29 @@ END MODULE MODI_STATION_n
! A. Lemonsu 19/11/2002
! P. Aumond 01/07/2011: add model levels
! C. Lac 04/2013: correction on the vertical levels
-! C. Lac 04/2013: add I/J positioning
+! C. Lac 04/2013: add I/JK positioning
! P. Wautelet 28/03/2018: replace TEMPORAL_DIST by DATETIME_DISTANCE
! P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
! P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
! R. Schoetter 11/2019: use LCARTESIAN instead of LSTATLAT for multiproc in cartesian
-! P. Wautelet 07/04/2022: rewrite types for stations
+! P. Wautelet 04/2022: restructure stations for better performance, reduce memory usage and correct some problems/bugs
!
! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
-USE MODD_CONF
-USE MODD_CST
+USE MODD_CONF, ONLY: LCARTESIAN
+USE MODD_CST, ONLY: XPI
USE MODD_DIAG_IN_RUN
-USE MODD_GRID
-USE MODD_PARAMETERS
+USE MODD_GRID, ONLY: XBETA, XLON0, XRPK
+USE MODD_PARAMETERS, ONLY: JPVEXT, XUNDEF
USE MODD_PARAM_n, ONLY: CRAD
USE MODD_STATION_n
USE MODD_ALLSTATION_n, ONLY: LDIAG_SURFRAD
-USE MODD_TIME, ONLY: tdtexp
USE MODD_TIME_n, ONLY: tdtcur
!
-USE MODE_ll
-!
-USE MODI_WATER_SUM
-!
+USE MODE_STATION_TOOLS, ONLY: STATION_INTERP_2D, STATION_INTERP_2D_U, STATION_INTERP_2D_V
!
!
IMPLICIT NONE
@@ -129,21 +125,6 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PP ! pressure
! 0.2 declaration of local variables
!
!
-INTEGER :: IIB ! current processor domain sizes
-INTEGER :: IJB !
-INTEGER :: IIE !
-INTEGER :: IJE !
-INTEGER :: IIU !
-INTEGER :: IJU !
-!
-REAL, DIMENSION(SIZE(PXHAT)) :: ZXHATM ! mass point coordinates
-REAL, DIMENSION(SIZE(PYHAT)) :: ZYHATM ! mass point coordinates
-!
-REAL, DIMENSION(SIZE(PSV,1),SIZE(PSV,2),SIZE(PSV,3),SIZE(PSV,4)) :: ZWORK !
-!
-LOGICAL :: GSTORE ! storage occurs at this time step
-!
-!
INTEGER :: IN ! time index
INTEGER :: JSV ! loop counter
!
@@ -151,33 +132,8 @@ 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 :: IRESP ! return code
-INTEGER :: I ! loop for stations
-INTEGER :: J ! loop for levels
-
-!
-!----------------------------------------------------------------------------
-!
-!* 2. PRELIMINARIES
-! -------------
-!
-!* 2.1 Indices
-! -------
-!
-CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
-!
-!
-!* 2.2 Interpolations of model variables to mass points
-! ------------------------------------------------
-!
-IIU=SIZE(PXHAT)
-IJU=SIZE(PYHAT)
-!
-ZXHATM(1:IIU-1)=0.5*PXHAT(1:IIU-1)+0.5*PXHAT(2:IIU )
-ZXHATM( IIU )=1.5*PXHAT( IIU )-0.5*PXHAT( IIU-1)
-!
-ZYHATM(1:IJU-1)=0.5*PYHAT(1:IJU-1)+0.5*PYHAT(2:IJU )
-ZYHATM( IJU )=1.5*PYHAT( IJU )-0.5*PYHAT( IJU-1)
+INTEGER :: JS ! loop for stations
+INTEGER :: JK ! loop for levels
!
!----------------------------------------------------------------------------
!
@@ -189,260 +145,73 @@ IF ( TSTATIONS_TIME%XTIME_CUR == XUNDEF ) TSTATIONS_TIME%XTIME_CUR = TSTATIONS_T
TSTATIONS_TIME%XTIME_CUR = TSTATIONS_TIME%XTIME_CUR + PTSTEP
!
IF ( TSTATIONS_TIME%XTIME_CUR >= TSTATIONS_TIME%XTSTEP - 1.E-10 ) THEN
- GSTORE = .TRUE.
- TSTATIONS_TIME%XTIME_CUR = TSTATIONS_TIME%XTIME_CUR - TSTATIONS_TIME%XTSTEP
- TSTATIONS_TIME%N_CUR = TSTATIONS_TIME%N_CUR + 1
- IN = TSTATIONS_TIME%N_CUR
-ELSE
- GSTORE = .FALSE.
-END IF
-!
-IF (GSTORE) THEN
-#if 0
- tstations_time%tpdates(in)%date%year = tdtexp%date%year
- tstations_time%tpdates(in)%date%month = tdtexp%date%month
- tstations_time%tpdates(in)%date%day = tdtexp%date%day
- tstations_time%tpdates(in)%xtime = tdtexp%xtime + ( in - 1 ) * tstation%step
-#else
+ TSTATIONS_TIME%XTIME_CUR = TSTATIONS_TIME%XTIME_CUR - TSTATIONS_TIME%XTSTEP
+ TSTATIONS_TIME%N_CUR = TSTATIONS_TIME%N_CUR + 1
+ IN = TSTATIONS_TIME%N_CUR
tstations_time%tpdates(in) = tdtcur
-#endif
+ELSE
+ !No station storage at this time step
+ RETURN
END IF
!
-!
!----------------------------------------------------------------------------
!
!* 8. DATA RECORDING
! --------------
!
-IF (GSTORE) THEN
- DO I=1,NUMBSTAT
- !
- 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
- J=1
- DO WHILE ((STATION_INTERP_2D(PZ(:,:,J))-STATION_INTERP_2D(PZ(:,:,2))) &
- < TSTATIONS(I)%XZ)
- J = J + 1
- END DO
- IF (((STATION_INTERP_2D(PZ(:,:,J))-STATION_INTERP_2D(PZ(:,:,2)))-TSTATIONS(I)%XZ)>&
- (TSTATIONS(I)%XZ-(STATION_INTERP_2D(PZ(:,:,J-1))-STATION_INTERP_2D(PZ(:,:,2))))) THEN
- J=J-1
- ENDIF
- END IF
- !
- IF (LCARTESIAN) THEN
- TSTATIONS(I)%XZON (IN) = STATION_INTERP_2D_U(PU(:,:,J))
- TSTATIONS(I)%XMER (IN) = STATION_INTERP_2D_V(PV(:,:,J))
- ELSE
- ZU_STAT = STATION_INTERP_2D_U(PU(:,:,J))
- ZV_STAT = STATION_INTERP_2D_V(PV(:,:,J))
- ZGAM = (XRPK * (TSTATIONS(I)%XLON - XLON0) - XBETA)*(XPI/180.)
- TSTATIONS(I)%XZON (IN) = ZU_STAT * COS(ZGAM) + ZV_STAT * SIN(ZGAM)
- TSTATIONS(I)%XMER (IN) = - ZU_STAT * SIN(ZGAM) + ZV_STAT * COS(ZGAM)
- ENDIF
- TSTATIONS(I)%XW (IN) = STATION_INTERP_2D(PW(:,:,J))
- TSTATIONS(I)%XTH (IN) = STATION_INTERP_2D(PTH(:,:,J))
- TSTATIONS(I)%XP (IN) = STATION_INTERP_2D(PP(:,:,J))
- !
- DO JSV=1,SIZE(PR,4)
- TSTATIONS(I)%XR (IN,JSV) = STATION_INTERP_2D(PR(:,:,J,JSV))
- END DO
- !
- DO JSV=1,SIZE(PSV,4)
- TSTATIONS(I)%XSV (IN,JSV) = STATION_INTERP_2D(PSV(:,:,J,JSV))
- END DO
- !
- IF (SIZE(PTKE)>0) TSTATIONS(I)%XTKE (IN) = STATION_INTERP_2D(PTKE(:,:,J))
- IF (SIZE(PTS) >0) TSTATIONS(I)%XTSRAD(IN) = STATION_INTERP_2D(PTS)
- TSTATIONS(I)%XZS = STATION_INTERP_2D(PZ(:,:,1+JPVEXT))
- !
- IF (LDIAG_SURFRAD) THEN
- TSTATIONS(I)%XZON10M(IN) = STATION_INTERP_2D(XCURRENT_ZON10M)
- TSTATIONS(I)%XMER10M(IN) = STATION_INTERP_2D(XCURRENT_MER10M)
- TSTATIONS(I)%XT2M (IN) = STATION_INTERP_2D(XCURRENT_T2M )
- TSTATIONS(I)%XQ2M (IN) = STATION_INTERP_2D(XCURRENT_Q2M )
- TSTATIONS(I)%XHU2M (IN) = STATION_INTERP_2D(XCURRENT_HU2M )
- TSTATIONS(I)%XRN (IN) = STATION_INTERP_2D(XCURRENT_RN )
- TSTATIONS(I)%XH (IN) = STATION_INTERP_2D(XCURRENT_H )
- TSTATIONS(I)%XLE (IN) = STATION_INTERP_2D(XCURRENT_LE )
- TSTATIONS(I)%XLEI (IN) = STATION_INTERP_2D(XCURRENT_LEI )
- TSTATIONS(I)%XGFLUX (IN) = STATION_INTERP_2D(XCURRENT_GFLUX )
- IF (CRAD /= 'NONE') THEN
- TSTATIONS(I)%XSWD (IN) = STATION_INTERP_2D(XCURRENT_SWD )
- TSTATIONS(I)%XSWU (IN) = STATION_INTERP_2D(XCURRENT_SWU )
- TSTATIONS(I)%XLWD (IN) = STATION_INTERP_2D(XCURRENT_LWD )
- TSTATIONS(I)%XLWU (IN) = STATION_INTERP_2D(XCURRENT_LWU )
- TSTATIONS(I)%XSWDIR (IN) = STATION_INTERP_2D(XCURRENT_SWDIR )
- TSTATIONS(I)%XSWDIFF(IN) = STATION_INTERP_2D(XCURRENT_SWDIFF)
- TSTATIONS(I)%XDSTAOD(IN) = STATION_INTERP_2D(XCURRENT_DSTAOD)
- ENDIF
- TSTATIONS(I)%XSFCO2 (IN) = STATION_INTERP_2D(XCURRENT_SFCO2 )
- ENDIF
+STATION: DO JS = 1,NUMBSTAT_LOC
+ JK = TSTATIONS(JS)%NK
+
+ IF (LCARTESIAN) THEN
+ TSTATIONS(JS)%XZON(IN) = STATION_INTERP_2D_U( TSTATIONS(JS), PU(:,:,JK) )
+ TSTATIONS(JS)%XMER(IN) = STATION_INTERP_2D_V( TSTATIONS(JS), PV(:,:,JK) )
+ ELSE
+ ZU_STAT = STATION_INTERP_2D_U( TSTATIONS(JS), PU(:,:,JK) )
+ ZV_STAT = STATION_INTERP_2D_V( TSTATIONS(JS), PV(:,:,JK) )
+ ZGAM = (XRPK * (TSTATIONS(JS)%XLON - XLON0) - XBETA)*(XPI/180.)
+ TSTATIONS(JS)%XZON(IN) = ZU_STAT * COS(ZGAM) + ZV_STAT * SIN(ZGAM)
+ TSTATIONS(JS)%XMER(IN) = - ZU_STAT * SIN(ZGAM) + ZV_STAT * COS(ZGAM)
+ END IF
+ TSTATIONS(JS)%XW (IN) = STATION_INTERP_2D( TSTATIONS(JS), PW(:,:,JK) )
+ TSTATIONS(JS)%XTH(IN) = STATION_INTERP_2D( TSTATIONS(JS), PTH(:,:,JK) )
+ TSTATIONS(JS)%XP (IN) = STATION_INTERP_2D( TSTATIONS(JS), PP(:,:,JK) )
- !
- END IF
-!
-!----------------------------------------------------------------------------
-!
-!* 11. EXCHANGE OF INFORMATION BETWEEN PROCESSORS
-! ------------------------------------------
-!
-!* 11.2 data stored
-! -----------
-!
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XX )
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XY )
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XZ )
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XLON )
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XLAT )
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XZON (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XMER (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XW (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XP (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XTH (IN))
DO JSV=1,SIZE(PR,4)
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XR (IN,JSV))
+ TSTATIONS(JS)%XR(IN,JSV) = STATION_INTERP_2D( TSTATIONS(JS), PR(:,:,JK,JSV) )
END DO
+
DO JSV=1,SIZE(PSV,4)
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XSV (IN,JSV))
+ TSTATIONS(JS)%XSV(IN,JSV) = STATION_INTERP_2D( TSTATIONS(JS), PSV(:,:,JK,JSV) )
END DO
- IF (SIZE(PTKE)>0) CALL DISTRIBUTE_STATION(TSTATIONS(I)%XTKE (IN))
- IF (SIZE(PTS) >0) CALL DISTRIBUTE_STATION(TSTATIONS(I)%XTSRAD(IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XZS )
- IF (LDIAG_SURFRAD) THEN
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XT2M (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XQ2M (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XHU2M (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XZON10M (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XMER10M (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XRN (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XH (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XLE (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XLEI (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XGFLUX (IN))
- IF (CRAD /= 'NONE') THEN
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XSWD (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XSWU (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XLWD (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XLWU (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XSWDIR (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XSWDIFF (IN))
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XDSTAOD (IN))
- END IF
- CALL DISTRIBUTE_STATION(TSTATIONS(I)%XSFCO2 (IN))
- ENDIF
- ENDDO
- !
-END IF
-!
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-!
-CONTAINS
-!
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-!
-FUNCTION STATION_INTERP_2D(PA) RESULT(PB)
-!
-REAL, DIMENSION(:,:), INTENT(IN) :: PA
-REAL :: PB
-!
-INTEGER :: JI, JJ
-!
-IF (SIZE(PA,1)==2) THEN
- JI=1
- JJ=1
-ELSEIF (L1D) THEN
- JI=2
- JJ=2
-ELSE
- 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.-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
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-! MODIFS
-FUNCTION STATION_INTERP_2D_U(PA) RESULT(PB)
-!
-REAL, DIMENSION(:,:), INTENT(IN) :: PA
-REAL :: PB
-!
-INTEGER :: JI, JJ
-!
-IF (SIZE(PA,1)==2) THEN
- JI=1
- JJ=1
-ELSEIF (L1D) THEN
- JI=2
- JJ=2
-ELSE
- 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.- 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
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-! MODIFS
-FUNCTION STATION_INTERP_2D_V(PA) RESULT(PB)
-!
-REAL, DIMENSION(:,:), INTENT(IN) :: PA
-REAL :: PB
-!
-INTEGER :: JI, JJ
-!
-IF (SIZE(PA,1)==2) THEN
- JI=1
- JJ=1
-ELSEIF (L1D) THEN
- JI=2
- JJ=2
-ELSE
- 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.- 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
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-SUBROUTINE DISTRIBUTE_STATION(PAS)
-!
-REAL, INTENT(INOUT) :: PAS
-!
-INTEGER :: IINFO_ll ! return code
+ IF (SIZE(PTKE)>0) TSTATIONS(JS)%XTKE(IN) = STATION_INTERP_2D( TSTATIONS(JS), PTKE(:,:,JK) )
+ IF ( CRAD /= 'NONE' ) TSTATIONS(JS)%XTSRAD(IN) = STATION_INTERP_2D( TSTATIONS(JS), PTS )
+ TSTATIONS(JS)%XZS = STATION_INTERP_2D( TSTATIONS(JS), PZ(:,:,1+JPVEXT))
-IF ( .NOT. TSTATIONS(I)%LPRESENT ) PAS = 0.
-CALL REDUCESUM_ll(PAS,IINFO_ll)
+ IF ( LDIAG_SURFRAD ) THEN
+ TSTATIONS(JS)%XZON10M(IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_ZON10M )
+ TSTATIONS(JS)%XMER10M(IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_MER10M )
+ TSTATIONS(JS)%XT2M (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_T2M )
+ TSTATIONS(JS)%XQ2M (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_Q2M )
+ TSTATIONS(JS)%XHU2M (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_HU2M )
+ TSTATIONS(JS)%XRN (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_RN )
+ TSTATIONS(JS)%XH (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_H )
+ TSTATIONS(JS)%XLE (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_LE )
+ TSTATIONS(JS)%XLEI (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_LEI )
+ TSTATIONS(JS)%XGFLUX (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_GFLUX )
+ IF ( CRAD /= 'NONE' ) THEN
+ TSTATIONS(JS)%XSWD (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_SWD )
+ TSTATIONS(JS)%XSWU (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_SWU )
+ TSTATIONS(JS)%XLWD (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_LWD )
+ TSTATIONS(JS)%XLWU (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_LWU )
+ TSTATIONS(JS)%XSWDIR (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_SWDIR )
+ TSTATIONS(JS)%XSWDIFF(IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_SWDIFF )
+ TSTATIONS(JS)%XDSTAOD(IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_DSTAOD )
+ END IF
+ TSTATIONS(JS)%XSFCO2(IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_SFCO2 )
+ END IF
+END DO STATION
!
-END SUBROUTINE DISTRIBUTE_STATION
!----------------------------------------------------------------------------
!
END SUBROUTINE STATION_n
diff --git a/src/MNH/write_stationn.f90 b/src/MNH/write_stationn.f90
index 0746c64d2..9f65f84e1 100644
--- a/src/MNH/write_stationn.f90
+++ b/src/MNH/write_stationn.f90
@@ -61,7 +61,7 @@ END MODULE MODI_WRITE_STATION_n
! P. Wautelet 09/10/2020: Write_diachro: use new datatype tpfields
! P. Wautelet 03/03/2021: budgets: add tbudiachrometadata type (useful to pass more information to Write_diachro)
! P. Wautelet 04/02/2022: use TSVLIST to manage metadata of scalar variables
-! P. Wautelet 07/04/2022: rewrite types for stations
+! P. Wautelet 04/2022: restructure stations for better performance, reduce memory usage and correct some problems/bugs
! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -70,17 +70,22 @@ END MODULE MODI_WRITE_STATION_n
USE MODD_ALLSTATION_n, ONLY: LDIAG_SURFRAD
use MODD_BUDGET, ONLY: tbudiachrometadata
USE MODD_CONF, ONLY: LCARTESIAN
+USE MODD_CONF_n, ONLY: NRR
USE MODD_CST, ONLY: XRV
-USE MODD_IO, ONLY: TFILEDATA
+USE MODD_IO, ONLY: ISNPROC, ISP, TFILEDATA
+USE MODD_MPIF
USE MODD_NSV, ONLY: tsvlist, nsv, nsv_aer, nsv_aerbeg, nsv_aerend, &
nsv_dst, nsv_dstbeg, nsv_dstend, nsv_slt, nsv_sltbeg, nsv_sltend
-USE MODD_PARAM_n, ONLY: CRAD, CSURF
+USE MODD_PARAM_n, ONLY: CRAD, CSURF, CTURB
USE MODD_PARAMETERS, ONLY: XUNDEF
-USE MODD_STATION_n, only: NUMBSTAT, TSTATIONS
+USE MODD_PRECISION, ONLY: MNHINT_MPI, MNHREAL_MPI
+USE MODD_STATION_n, only: NUMBSTAT_LOC, TSTATIONS, tstations_time
+USE MODD_TYPE_STATION, ONLY: TSTATIONDATA
!
USE MODE_AERO_PSD
USE MODE_DUST_PSD
USE MODE_SALT_PSD
+USE MODE_STATION_TOOLS, ONLY: STATION_ALLOCATE
use MODE_WRITE_DIACHRO, ONLY: Write_diachro
!
IMPLICIT NONE
@@ -94,13 +99,199 @@ TYPE(TFILEDATA), INTENT(IN) :: TPDIAFILE ! diachronic file to write
!
! 0.2 declaration of local variables
!
-INTEGER :: II ! loop
+INTEGER, PARAMETER :: ITAG = 100
+INTEGER :: IERR
+INTEGER :: JP, JS
+INTEGER :: IDX
+INTEGER :: INUMSTAT ! Total number of stations (for the current model)
+INTEGER :: IPACKSIZE ! Size of the ZPACK buffer
+INTEGER :: IPOS ! Position in the ZPACK buffer
+INTEGER :: ISTORE
+INTEGER, DIMENSION(:), ALLOCATABLE :: INSTATPRC ! Array to store the number of stations per process (for the current model)
+INTEGER, DIMENSION(:), ALLOCATABLE :: ISTATIDS ! Intermediate array for MPI communication
+INTEGER, DIMENSION(:), ALLOCATABLE :: ISTATPRCRANK ! Array to store the ranks of the processes where the stations are
+INTEGER, DIMENSION(:), ALLOCATABLE :: IDS ! Array to store the station number to send
+INTEGER, DIMENSION(:), ALLOCATABLE :: IDISP ! Array to store the displacements for MPI communications
+REAL, DIMENSION(:), ALLOCATABLE :: ZPACK ! Buffer to store raw data of a station (used for MPI communication)
+TYPE(TSTATIONDATA) :: TZSTATION
!
!----------------------------------------------------------------------------
-!
-DO II = 1, NUMBSTAT
- CALL STATION_DIACHRO_n( TSTATIONS(II) )
-ENDDO
+
+ALLOCATE( INSTATPRC(ISNPROC) )
+ALLOCATE( IDS(NUMBSTAT_LOC) )
+
+!Gather number of station present on each process
+CALL MPI_ALLGATHER( NUMBSTAT_LOC, 1, MNHINT_MPI, INSTATPRC, 1, MNHINT_MPI, TPDIAFILE%NMPICOMM, IERR )
+
+!Store the identification number of local stations (these numbers are globals)
+DO JS = 1, NUMBSTAT_LOC
+ IDS(JS) = TSTATIONS(JS)%NID
+END DO
+
+ALLOCATE( IDISP(ISNPROC) )
+IDISP(1) = 0
+DO JP = 2, ISNPROC
+ IDISP(JP) = IDISP(JP-1) + INSTATPRC(JP-1)
+END DO
+
+INUMSTAT = SUM( INSTATPRC(:) )
+ALLOCATE( ISTATIDS(INUMSTAT) )
+ALLOCATE( ISTATPRCRANK(INUMSTAT) )
+
+!Gather the list of all the stations of all processes
+CALL MPI_ALLGATHERV( IDS(:), NUMBSTAT_LOC, MNHINT_MPI, ISTATIDS(:), INSTATPRC(:), &
+ IDISP(:), MNHINT_MPI, TPDIAFILE%NMPICOMM, IERR )
+
+!Store the rank of each process corresponding to a given station
+IDX = 1
+ISTATPRCRANK(:) = -1
+DO JP = 1, ISNPROC
+ DO JS = 1, INSTATPRC(JP)
+ ISTATPRCRANK(ISTATIDS(IDX)) = JP
+ IDX = IDX + 1
+ END DO
+END DO
+
+CALL STATION_ALLOCATE( TZSTATION, SIZE( tstations_time%tpdates ) )
+
+!Determine the size of the ZPACK buffer used to transfer station data in 1 MPI communication
+IF ( ISNPROC > 1 ) THEN
+ ISTORE = SIZE( TSTATIONS_TIME%TPDATES )
+ IPACKSIZE = 7
+ IPACKSIZE = IPACKSIZE + ISTORE * ( 5 + NRR + NSV )
+ IF ( CTURB == 'TKEL') IPACKSIZE = IPACKSIZE + ISTORE !Tke term
+ IF ( CRAD /= 'NONE' ) IPACKSIZE = IPACKSIZE + ISTORE !XTSRAD term
+ IF (LDIAG_SURFRAD) THEN
+ IF ( CSURF == 'EXTE' ) IPACKSIZE = IPACKSIZE + ISTORE * 10
+ IF ( CRAD /= 'NONE' ) IPACKSIZE = IPACKSIZE + ISTORE * 7
+ IPACKSIZE = IPACKSIZE + ISTORE !XSFCO2 term
+ END IF
+
+ ALLOCATE( ZPACK(IPACKSIZE) )
+END IF
+
+IDX = 1
+
+STATION: DO JS = 1, INUMSTAT
+ IF ( ISTATPRCRANK(JS) == TPDIAFILE%NMASTER_RANK ) THEN
+ !No communication necessary, the station data is already on the writer process
+ IF ( ISP == TPDIAFILE%NMASTER_RANK ) THEN
+ TZSTATION = TSTATIONS(IDX)
+ IDX = IDX + 1
+ END IF
+ ELSE
+ !The station data is not on the writer process
+ IF ( ISP == ISTATPRCRANK(JS) ) THEN
+ ! This process has the data and needs to send it to the writer process
+ IPOS = 1
+ ZPACK(IPOS) = TSTATIONS(IDX)%NID; IPOS = IPOS + 1
+ ZPACK(IPOS) = TSTATIONS(IDX)%XX; IPOS = IPOS + 1
+ ZPACK(IPOS) = TSTATIONS(IDX)%XY; IPOS = IPOS + 1
+ ZPACK(IPOS) = TSTATIONS(IDX)%XZ; IPOS = IPOS + 1
+ ZPACK(IPOS) = TSTATIONS(IDX)%XLON; IPOS = IPOS + 1
+ ZPACK(IPOS) = TSTATIONS(IDX)%XLAT; IPOS = IPOS + 1
+ ZPACK(IPOS) = TSTATIONS(IDX)%XZS; IPOS = IPOS + 1
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XZON(:); IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XMER(:); IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XW(:); IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XP(:); IPOS = IPOS + ISTORE
+ IF ( CTURB == 'TKEL') THEN
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XTKE(:); IPOS = IPOS + ISTORE
+ END IF
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XTH(:); IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE*NRR-1) = RESHAPE( TSTATIONS(IDX)%XR(:,:), [ISTORE*NRR] ); IPOS = IPOS + ISTORE * NRR
+ ZPACK(IPOS:IPOS+ISTORE*NSV-1) = RESHAPE( TSTATIONS(IDX)%XSV(:,:), [ISTORE*NSV] ); IPOS = IPOS + ISTORE * NSV
+ IF ( CRAD /= 'NONE' ) THEN
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XTSRAD(:); IPOS = IPOS + ISTORE
+ END IF
+ IF (LDIAG_SURFRAD) THEN
+ IF ( CSURF == 'EXTE') THEN
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XT2M; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XQ2M; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XHU2M; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XZON10M; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XMER10M; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XRN; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XH; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XLE; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XGFLUX; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XLEI; IPOS = IPOS + ISTORE
+ END IF
+ IF ( CRAD /= 'NONE' ) THEN
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XSWD; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XSWU; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XLWD; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XLWU; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XSWDIR; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XSWDIFF; IPOS = IPOS + ISTORE
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XDSTAOD; IPOS = IPOS + ISTORE
+ END IF
+ ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XSFCO2; IPOS = IPOS + ISTORE
+ END IF
+
+ CALL MPI_SEND( TSTATIONS(IDX)%CNAME, LEN(TSTATIONS(IDX)%CNAME), MPI_CHARACTER, TPDIAFILE%NMASTER_RANK - 1, &
+ ITAG, TPDIAFILE%NMPICOMM, IERR )
+ CALL MPI_SEND( ZPACK, IPACKSIZE, MNHREAL_MPI, TPDIAFILE%NMASTER_RANK - 1, ITAG, TPDIAFILE%NMPICOMM, IERR )
+
+ IDX = IDX + 1
+
+ ELSE IF ( ISP == TPDIAFILE%NMASTER_RANK ) THEN
+ ! This process is the writer and will receive the station data from its owner
+ CALL MPI_RECV( TZSTATION%CNAME, LEN(TZSTATION%CNAME), MPI_CHARACTER, &
+ ISTATPRCRANK(JS) - 1, ITAG, TPDIAFILE%NMPICOMM, MPI_STATUS_IGNORE, IERR )
+ CALL MPI_RECV( ZPACK, IPACKSIZE, MNHREAL_MPI, ISTATPRCRANK(JS) - 1, ITAG, TPDIAFILE%NMPICOMM, MPI_STATUS_IGNORE, IERR )
+
+ IPOS = 1
+ TZSTATION%NID = NINT( ZPACK(IPOS) ); IPOS = IPOS + 1
+ TZSTATION%XX = ZPACK(IPOS); IPOS = IPOS + 1
+ TZSTATION%XY = ZPACK(IPOS); IPOS = IPOS + 1
+ TZSTATION%XZ = ZPACK(IPOS); IPOS = IPOS + 1
+ TZSTATION%XLON = ZPACK(IPOS); IPOS = IPOS + 1
+ TZSTATION%XLAT = ZPACK(IPOS); IPOS = IPOS + 1
+ TZSTATION%XZS = ZPACK(IPOS); IPOS = IPOS + 1
+ TZSTATION%XZON(:) = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XMER(:) = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XW(:) = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XP(:) = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ IF ( CTURB == 'TKEL') THEN
+ TZSTATION%XTKE(:) = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ END IF
+ TZSTATION%XTH(:) = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XR(:,:) = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*NRR-1), [ ISTORE, NRR ] ); IPOS = IPOS + ISTORE * NRR
+ TZSTATION%XSV(:,:) = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*NSV-1), [ ISTORE, NSV ] ); IPOS = IPOS + ISTORE * NSV
+ IF ( CRAD /= 'NONE' ) THEN
+ TZSTATION%XTSRAD(:) = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ END IF
+ IF (LDIAG_SURFRAD) THEN
+ IF ( CSURF == 'EXTE' ) THEN
+ TZSTATION%XT2M = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XQ2M = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XHU2M = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XZON10M = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XMER10M = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XRN = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XH = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XLE = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XGFLUX = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XLEI = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ END IF
+ IF ( CRAD /= 'NONE' ) THEN
+ TZSTATION%XSWD = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XSWU = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XLWD = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XLWU = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XSWDIR = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XSWDIFF = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ TZSTATION%XDSTAOD = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ END IF
+ TZSTATION%XSFCO2 = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
+ END IF
+ END IF
+ END IF
+
+ CALL STATION_DIACHRO_n( TZSTATION )
+
+END DO STATION
!
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
@@ -145,13 +336,10 @@ type(tbudiachrometadata) :: tzbudiachro
type(tfieldmetadata_base), dimension(:), allocatable :: tzfields
!
!----------------------------------------------------------------------------
-IF (TPSTATION%XX==XUNDEF) RETURN
-IF (TPSTATION%XY==XUNDEF) RETURN
!
IPROC = 8 + SIZE(TPSTATION%XR,2) + SIZE(TPSTATION%XSV,2)
-IF (TPSTATION%XX==XUNDEF) IPROC = IPROC + 2
-IF (SIZE(TPSTATION%XTKE )>0) IPROC = IPROC + 1
+IF ( CTURB == 'TKEL' ) IPROC = IPROC + 1
IF (LDIAG_SURFRAD) THEN
IF(CSURF=="EXTE") IPROC = IPROC + 10
IF(CRAD/="NONE") IPROC = IPROC + 7
@@ -159,8 +347,8 @@ END IF
IF (LORILAM) IPROC = IPROC + JPMODE*(3+NSOA+NCARB+NSP)
IF (LDUST) IPROC = IPROC + NMODE_DST*3
IF (LSALT) IPROC = IPROC + NMODE_SLT*3
-IF (ANY(TPSTATION%XTSRAD(:)/=XUNDEF)) IPROC = IPROC + 1
-IF (ANY(TPSTATION%XSFCO2(:)/=XUNDEF)) IPROC = IPROC + 1
+IF ( CRAD /= 'NONE' ) IPROC = IPROC + 1
+IPROC = IPROC + 1 ! XSFCO2 term
!
ALLOCATE (ZWORK6(1,1,1,SIZE(tstations_time%tpdates),1,IPROC))
ALLOCATE (YCOMMENT(IPROC))
@@ -390,7 +578,7 @@ DO JRR=1,SIZE(TPSTATION%XR,2)
END IF
END DO
!
-IF (SIZE(TPSTATION%XTKE,1)>0) THEN
+IF ( CTURB == 'TKEL' ) THEN
JPROC = JPROC+1
YTITLE (JPROC) = 'Tke'
YUNIT (JPROC) = 'm2 s-2'
@@ -642,7 +830,7 @@ IF (SIZE(TPSTATION%XSV,2)>=1) THEN
END IF
END IF
-IF (ANY(TPSTATION%XTSRAD(:)/=XUNDEF)) THEN
+IF ( CRAD /= 'NONE' ) THEN
JPROC = JPROC+1
YTITLE (JPROC) = 'Tsrad'
YUNIT (JPROC) = 'K'
@@ -650,13 +838,13 @@ IF (ANY(TPSTATION%XTSRAD(:)/=XUNDEF)) THEN
ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XTSRAD(:)
END IF
!
-IF (ANY(TPSTATION%XSFCO2(:)/=XUNDEF)) THEN
+! IF (ANY(TPSTATION%XSFCO2(:)/=XUNDEF)) THEN
JPROC = JPROC+1
YTITLE (JPROC) = 'SFCO2'
YUNIT (JPROC) = 'mg m-2 s-1'
YCOMMENT (JPROC) = 'CO2 Surface Flux'
ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XSFCO2(:)
-END IF
+! END IF
!
!----------------------------------------------------------------------------
!
@@ -739,7 +927,5 @@ DEALLOCATE (YUNIT )
DEALLOCATE (IGRID )
!----------------------------------------------------------------------------
END SUBROUTINE STATION_DIACHRO_n
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
!
END SUBROUTINE WRITE_STATION_n
--
GitLab