From 446ef213fc9a3781f8cb70af609705884ff2d0b9 Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Tue, 22 Mar 2022 13:52:16 +0100
Subject: [PATCH] Philippe 22/03/2022: LES budgets: LES averaging periods are
more reliable (compute with integers instead of reals)
(cherry picked from commit 1ca4e62cf481d02e29ecf61f39173e05d5a70cec)
---
src/MNH/ini_lesn.f90 | 41 ++++++++++++++++++++++++++
src/MNH/modd_lesn.f90 | 15 +++++++++-
src/MNH/mode_les_diachro.f90 | 57 ++++++++++--------------------------
3 files changed, 70 insertions(+), 43 deletions(-)
diff --git a/src/MNH/ini_lesn.f90 b/src/MNH/ini_lesn.f90
index 05a30c423..84f90dba5 100644
--- a/src/MNH/ini_lesn.f90
+++ b/src/MNH/ini_lesn.f90
@@ -40,6 +40,7 @@
! P. Wautelet 04/01/2021: bugfix: nles_k was used instead of nspectra_k for a loop index
! P. Wautelet 30/03/2021: budgets: LES cartesian subdomain limits are defined in the physical domain
! P. Wautelet 09/07/2021: bugfix: altitude levels are on the correct grid position (mass point)
+! P. Wautelet 22/03/2022: LES averaging periods are more reliable (compute with integers instead of reals)
! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -342,6 +343,46 @@ IF (NLES_TIMES==0) THEN
RETURN
END IF
!
+!* 3.8 Averaging
+! ---------
+IF ( XLES_TEMP_MEAN_END == XUNDEF &
+ .OR. XLES_TEMP_MEAN_START == XUNDEF &
+ .OR. XLES_TEMP_MEAN_STEP == XUNDEF ) THEN
+ !No LES temporal averaging
+ NLES_MEAN_TIMES = 0
+ NLES_MEAN_STEP = NNEGUNDEF
+ NLES_MEAN_START = NNEGUNDEF
+ NLES_MEAN_END = NNEGUNDEF
+ELSE
+ !LES temporal averaging is enabled
+ !Ensure that XLES_TEMP_MEAN_END is not after segment end
+ XLES_TEMP_MEAN_END = MIN( XLES_TEMP_MEAN_END, XSEGLEN - DYN_MODEL(1)%XTSTEP )
+
+ NLES_MEAN_START = NINT( XLES_TEMP_MEAN_START / XTSTEP )
+
+ IF ( MODULO( NLES_MEAN_START, NLES_DTCOUNT ) /= 0 ) THEN
+ CMNHMSG(1) = 'XLES_TEMP_MEAN_START is not a multiple of XLES_TEMP_SAMPLING'
+ CMNHMSG(2) = 'LES averaging periods could be wrong'
+ CALL Print_msg( NVERB_WARNING, 'IO', 'INI_LES_n' )
+ END IF
+
+ NLES_MEAN_END = NINT( XLES_TEMP_MEAN_END / XTSTEP )
+
+ NLES_MEAN_STEP = NINT( XLES_TEMP_MEAN_STEP / XTSTEP )
+
+ IF ( NLES_MEAN_STEP < NLES_DTCOUNT ) &
+ CALL Print_msg( NVERB_ERROR, 'IO', 'INI_LES_n', 'XLES_TEMP_MEAN_STEP < XLES_TEMP_SAMPLING not allowed' )
+
+ IF ( MODULO( NLES_MEAN_STEP, NLES_DTCOUNT ) /= 0 ) THEN
+ CMNHMSG(1) = 'XLES_TEMP_MEAN_STEP is not a multiple of XLES_TEMP_SAMPLING'
+ CMNHMSG(2) = 'LES averaging periods could be wrong'
+ CALL Print_msg( NVERB_WARNING, 'IO', 'INI_LES_n' )
+ END IF
+
+ NLES_MEAN_TIMES = ( NLES_MEAN_END - NLES_MEAN_START ) / NLES_MEAN_STEP
+ !Add 1 averaging period if the last one is incomplete (for example: start=0., end=10., step=3.)
+ IF ( MODULO( NLES_MEAN_END - NLES_MEAN_START, NLES_MEAN_STEP ) > 0 ) NLES_MEAN_TIMES = NLES_MEAN_TIMES + 1
+END IF
!-------------------------------------------------------------------------------
!
!* 4. Number of vertical levels for local diagnostics
diff --git a/src/MNH/modd_lesn.f90 b/src/MNH/modd_lesn.f90
index 28db43c4d..c7b54dd8d 100644
--- a/src/MNH/modd_lesn.f90
+++ b/src/MNH/modd_lesn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1995-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1995-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.
@@ -42,6 +42,7 @@
! P. Wautelet 08/02/2019: add missing NULL association for pointers
! C. Lac 02/2019: add rain fraction as a LES diagnostic
! P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
+! P. Wautelet 22/03/2022: LES averaging periods are more reliable (compute with integers instead of reals)
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -61,6 +62,10 @@ TYPE LES_t
INTEGER :: NLES_TIMES ! number of LES computations in time
INTEGER :: NLES_DTCOUNT ! number of time steps between two LES comp.
INTEGER :: NLES_TCOUNT ! current time counter for LES comp.
+ INTEGER :: NLES_MEAN_TIMES ! Number of LES averaging periods
+ INTEGER :: NLES_MEAN_STEP ! number of time steps between two LES average comp.
+ INTEGER :: NLES_MEAN_START ! First time step number taken into account for LES averaging
+ INTEGER :: NLES_MEAN_END ! Last time step number taken into account for LES averaging
INTEGER :: NSPECTRA_NI ! number of wave lengths in I direction
INTEGER :: NSPECTRA_NJ ! number of wave lengths in J direction
!
@@ -660,6 +665,10 @@ TYPE(LES_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: LES_MODEL
INTEGER, POINTER :: NLES_TIMES=>NULL()
INTEGER, POINTER :: NLES_DTCOUNT=>NULL()
INTEGER, POINTER :: NLES_TCOUNT=>NULL()
+INTEGER, POINTER :: NLES_MEAN_TIMES => NULL()
+INTEGER, POINTER :: NLES_MEAN_STEP => NULL()
+INTEGER, POINTER :: NLES_MEAN_START => NULL()
+INTEGER, POINTER :: NLES_MEAN_END => NULL()
INTEGER, POINTER :: NSPECTRA_NI=>NULL()
INTEGER, POINTER :: NSPECTRA_NJ=>NULL()
type(date_time), dimension(:), pointer :: tles_dates => null()
@@ -1509,6 +1518,10 @@ LES_MODEL(KFROM)%XLES_RADEFF=>XLES_RADEFF
NLES_TIMES=>LES_MODEL(KTO)%NLES_TIMES
NLES_DTCOUNT=>LES_MODEL(KTO)%NLES_DTCOUNT
NLES_TCOUNT=>LES_MODEL(KTO)%NLES_TCOUNT
+NLES_MEAN_TIMES => LES_MODEL(KTO)%NLES_MEAN_TIMES
+NLES_MEAN_STEP => LES_MODEL(KTO)%NLES_MEAN_STEP
+NLES_MEAN_START => LES_MODEL(KTO)%NLES_MEAN_START
+NLES_MEAN_END => LES_MODEL(KTO)%NLES_MEAN_END
NSPECTRA_NI=>LES_MODEL(KTO)%NSPECTRA_NI
NSPECTRA_NJ=>LES_MODEL(KTO)%NSPECTRA_NJ
tles_dates=>les_model(kto)%tles_dates
diff --git a/src/MNH/mode_les_diachro.f90 b/src/MNH/mode_les_diachro.f90
index f56b2ce0f..7d8d69872 100644
--- a/src/MNH/mode_les_diachro.f90
+++ b/src/MNH/mode_les_diachro.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.
@@ -12,6 +12,7 @@
! P. Wautelet 10/2020: restructure subroutines to use tfield_metadata_base type
! P. Wautelet 03/03/2021: budgets: add tbudiachrometadata type (useful to pass more information to Write_diachro)
! P. Wautelet 11/03/2021: budgets: remove ptrajx/y/z optional dummy arguments of Write_diachro
+! P. Wautelet 22/03/2022: LES averaging periods are more reliable (compute with integers instead of reals)
!-----------------------------------------------------------------
!#######################
MODULE MODE_LES_DIACHRO
@@ -524,8 +525,7 @@ END SUBROUTINE LES_TIME_AVG
subroutine Les_time_avg_4D( pwork4, tpdates, kresp )
!########################################################
-use modd_les, only: nles_current_times, xles_temp_mean_start, xles_temp_mean_end, xles_temp_mean_step
-use modd_parameters, only: XUNDEF
+use modd_les_n, only: nles_dtcount, nles_mean_start, nles_mean_end, nles_mean_step, nles_mean_times, nles_times
use modd_time, only: tdtseg
use modd_type_date, only: date_time
@@ -537,60 +537,33 @@ real, dimension(:,:,:,:), allocatable, intent(inout) :: pwork4 ! cont
type(date_time), dimension(:), allocatable, intent(inout) :: tpdates
integer, intent(out) :: kresp ! return code (0 is ok)
!------------------------------------------------------------------------------
-integer :: jt ! time counter
-integer :: itime ! nb of avg. points
-integer :: iavg ! nb of avg. periods
integer :: javg ! loop counter on avg. periods
integer :: jk ! vertical loop counter
integer :: jp ! process loop counter
integer :: jsv ! scalar loop counter
-integer :: jx ! first spatial or spectral coordinate loop counter
-integer :: jy ! second spatial or spectral coordinate loop counter
integer :: jtb, jte
-real :: zles_temp_mean_start ! initial and end times
-real :: zles_temp_mean_end ! of one avergaing preiod
real, dimension(:,:,:,:), allocatable :: zwork4 ! contains averaged physical field
!------------------------------------------------------------------------------
-if ( xles_temp_mean_end == XUNDEF &
- .or. xles_temp_mean_start == XUNDEF &
- .or. xles_temp_mean_step == XUNDEF ) then
- kresp = -1
- return
-end if
-
-iavg = Int( xles_temp_mean_end - 1.e-10 - xles_temp_mean_start ) / xles_temp_mean_step + 1
-if ( iavg <= 0 ) then
+if ( nles_mean_times == 0 ) then
kresp = -1
return
end if
Deallocate( tpdates )
-Allocate( tpdates(iavg) )
-Allocate( zwork4(Size( pwork4, 1 ), iavg, Size( pwork4, 3 ), Size( pwork4, 4 )) )
+Allocate( tpdates(nles_mean_times) )
+Allocate( zwork4(Size( pwork4, 1 ), nles_mean_times, Size( pwork4, 3 ), Size( pwork4, 4 )) )
zwork4(:, :, :, :) = 0.
-do javg = 1, iavg
- zles_temp_mean_start = xles_temp_mean_start + (javg - 1) * xles_temp_mean_step
- zles_temp_mean_end = Min( xles_temp_mean_end, zles_temp_mean_start + xles_temp_mean_step )
-
- jtb = -1
- jte = -2
- do jt = 1, nles_current_times
- if ( xles_times(jt) >= zles_temp_mean_start ) then
- jtb = jt
- exit
- end if
- end do
- do jt = jtb, nles_current_times
- if ( xles_times(jt) <= zles_temp_mean_end ) then
- jte = jt
- else
- exit
- end if
- end do
+do javg = 1, nles_mean_times
+ jtb = ( nles_mean_start + ( javg - 1 ) * nles_mean_step ) / nles_dtcount
+ jte = MIN( jtb + nles_mean_step / nles_dtcount, nles_mean_end / nles_dtcount, nles_times )
+ ! jtb could be 0 if nles_mean_start is smaller than the first LES measurement
+ ! For example, it occurs if xles_temp_mean_start is smaller than xles_temp_sampling (if xles_temp_mean_start=0.)
+ ! Do this correction only after computation of jte
+ if ( jtb < 1 ) jtb = 1
do jp = 1, Size( pwork4, 4 )
do jsv = 1, Size( pwork4, 3 )
@@ -603,12 +576,12 @@ do javg = 1, iavg
tpdates(javg)%nyear = tdtseg%nyear
tpdates(javg)%nmonth = tdtseg%nmonth
tpdates(javg)%nday = tdtseg%nday
- tpdates(javg)%xtime = tdtseg%xtime + ( zles_temp_mean_start + zles_temp_mean_end ) / 2.
+ tpdates(javg)%xtime = tdtseg%xtime + ( xles_times(jtb) + xles_times(jte) ) / 2.
call Datetime_correctdate( tpdates(javg) )
end do
Deallocate( pwork4 )
-Allocate( pwork4(Size( zwork4, 1 ), iavg, Size( zwork4, 3 ), Size( zwork4, 4 )) )
+Allocate( pwork4(Size( zwork4, 1 ), nles_mean_times, Size( zwork4, 3 ), Size( zwork4, 4 )) )
pwork4 = zwork4
--
GitLab