diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index d8dc56c4a3a12cfd1d12951c6996139fe08e6554..25711a16726c5d06e5b4418456078eae8042fc79 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -1139,6 +1139,7 @@ ENDIF
!
IF (KMI == 1) THEN
LMEAN_FIELD = .FALSE.
+ LCOV_FIELD = .FALSE.
ENDIF
!
!-------------------------------------------------------------------------------
diff --git a/src/MNH/ini_mean_field.f90 b/src/MNH/ini_mean_field.f90
index 32ddfbb8847adcfe2cbe1887368abee74b2370b9..1abffffd053f5d6ead06c6979420e037acc58249 100644
--- a/src/MNH/ini_mean_field.f90
+++ b/src/MNH/ini_mean_field.f90
@@ -50,6 +50,7 @@ END MODULE MODI_INI_MEAN_FIELD
!! 10/2016 (C.Lac) Add max values
!! 02/2021 (T.Nagel) add passive scalar (XSVT) and UW wind component
!! 05/2021 (PA.Joulin) add wind turbine variables
+!! 11/2022 (E. Jezequel) add covariances
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -79,13 +80,18 @@ IF (CTURB /= 'NONE') XTKEM_MEAN = 0.0
XPABSM_MEAN = 0.0
XSVT_MEAN = 0.0
!
-XU2_MEAN = 0.0
-XV2_MEAN = 0.0
-XW2_MEAN = 0.0
-XUW_MEAN = 0.0
-XTH2_MEAN = 0.0
-XTEMP2_MEAN = 0.0
-XPABS2_MEAN = 0.0
+XU2_M2 = 0.0
+XV2_M2 = 0.0
+XW2_M2 = 0.0
+IF (LCOV_FIELD) THEN
+ XUV_MEAN = 0.0
+ XUW_MEAN = 0.0
+ XVW_MEAN = 0.0
+ XWTH_MEAN = 0.0
+END IF
+XTH2_M2 = 0.0
+XTEMP2_M2 = 0.0
+XPABS2_M2 = 0.0
!
IMI = GET_CURRENT_MODEL_INDEX()
IF (LMAIN_EOL .AND. IMI==NMODEL_EOL) THEN
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index 3d9aa05443f72ccb3c8beb9d0530ce48c56ed393..5d7ae53bf1d45b1e3fc868ca208593acd9762605 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -777,13 +777,21 @@ IF (LMEAN_FIELD) THEN
END IF
ALLOCATE(XPABSM_MEAN(IIU,IJU,IKU)) ; XPABSM_MEAN = 0.0
!
- ALLOCATE(XU2_MEAN(IIU,IJU,IKU)) ; XU2_MEAN = 0.0
- ALLOCATE(XV2_MEAN(IIU,IJU,IKU)) ; XV2_MEAN = 0.0
- ALLOCATE(XW2_MEAN(IIU,IJU,IKU)) ; XW2_MEAN = 0.0
- ALLOCATE(XUW_MEAN(IIU,IJU,IKU)) ; XUW_MEAN = 0.0
- ALLOCATE(XTH2_MEAN(IIU,IJU,IKU)) ; XTH2_MEAN = 0.0
- ALLOCATE(XTEMP2_MEAN(IIU,IJU,IKU)) ; XTEMP2_MEAN = 0.0
- ALLOCATE(XPABS2_MEAN(IIU,IJU,IKU)) ; XPABS2_MEAN = 0.0
+ ALLOCATE(XU2_M2(IIU,IJU,IKU)) ; XU2_M2 = 0.0
+!
+ ALLOCATE(XU2_M2(IIU,IJU,IKU)) ; XU2_M2 = 0.0
+ ALLOCATE(XV2_M2(IIU,IJU,IKU)) ; XV2_M2 = 0.0
+ ALLOCATE(XW2_M2(IIU,IJU,IKU)) ; XW2_M2 = 0.0
+ ALLOCATE(XTH2_M2(IIU,IJU,IKU)) ; XTH2_M2 = 0.0
+ ALLOCATE(XTEMP2_M2(IIU,IJU,IKU)) ; XTEMP2_M2 = 0.0
+ ALLOCATE(XPABS2_M2(IIU,IJU,IKU)) ; XPABS2_M2 = 0.0
+!
+ IF (LCOV_FIELD) THEN
+ ALLOCATE(XUV_MEAN(IIU,IJU,IKU)) ; XUV_MEAN = 0.0
+ ALLOCATE(XUW_MEAN(IIU,IJU,IKU)) ; XUW_MEAN = 0.0
+ ALLOCATE(XVW_MEAN(IIU,IJU,IKU)) ; XVW_MEAN = 0.0
+ ALLOCATE(XWTH_MEAN(IIU,IJU,IKU)) ; XWTH_MEAN = 0.0
+ END IF
!
ALLOCATE(XUM_MAX(IIU,IJU,IKU)) ; XUM_MAX = -1.E20
ALLOCATE(XVM_MAX(IIU,IJU,IKU)) ; XVM_MAX = -1.E20
@@ -798,6 +806,7 @@ IF (LMEAN_FIELD) THEN
END IF
ALLOCATE(XPABSM_MAX(IIU,IJU,IKU)) ; XPABSM_MAX = 0.0
ELSE
+!
ALLOCATE(XUM_MEAN(0,0,0))
ALLOCATE(XVM_MEAN(0,0,0))
ALLOCATE(XWM_MEAN(0,0,0))
@@ -807,13 +816,19 @@ ELSE
ALLOCATE(XTKEM_MEAN(0,0,0))
ALLOCATE(XPABSM_MEAN(0,0,0))
!
- ALLOCATE(XU2_MEAN(0,0,0))
- ALLOCATE(XV2_MEAN(0,0,0))
- ALLOCATE(XW2_MEAN(0,0,0))
- ALLOCATE(XUW_MEAN(0,0,0))
- ALLOCATE(XTH2_MEAN(0,0,0))
- ALLOCATE(XTEMP2_MEAN(0,0,0))
- ALLOCATE(XPABS2_MEAN(0,0,0))
+ ALLOCATE(XU2_M2(0,0,0))
+ ALLOCATE(XV2_M2(0,0,0))
+ ALLOCATE(XW2_M2(0,0,0))
+ ALLOCATE(XTH2_M2(0,0,0))
+ ALLOCATE(XTEMP2_M2(0,0,0))
+ ALLOCATE(XPABS2_M2(0,0,0))
+!
+ IF (LCOV_FIELD) THEN
+ ALLOCATE(XUV_MEAN(0,0,0))
+ ALLOCATE(XUW_MEAN(0,0,0))
+ ALLOCATE(XVW_MEAN(0,0,0))
+ ALLOCATE(XWTH_MEAN(0,0,0))
+ END IF
!
ALLOCATE(XUM_MAX(0,0,0))
ALLOCATE(XVM_MAX(0,0,0))
diff --git a/src/MNH/mean_field.f90 b/src/MNH/mean_field.f90
index 048eb689dfa04f798b16e555bb6b220e896b3484..a94d92633b762071be438dd28d4dd9c6f4bda493 100644
--- a/src/MNH/mean_field.f90
+++ b/src/MNH/mean_field.f90
@@ -48,6 +48,7 @@ END MODULE MODI_MEAN_FIELD
!! Original 07/2009
!! (C.Lac) 09/2016 Max values
!! (PA.Joulin) 12/2020 Wind turbine variables
+!! (E. Jezequel) 11/2022 Welford algorithm and covariances
!!---------------------------------------------------------------
!
!
@@ -67,6 +68,7 @@ USE MODD_EOL_SHARED_IO, ONLY: XTHRU_SUM, XTORQ_SUM, XPOW_SUM
USE MODD_EOL_ALM
USE MODD_EOL_ADNR
USE MODE_MODELN_HANDLER
+USE MODI_UPDATE_WELFORD
!
IMPLICIT NONE
@@ -82,6 +84,9 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PSVT
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZTEMPT
INTEGER :: IIU,IJU,IKU,IIB,IJB,IKB,IIE,IJE,IKE ! Arrays bounds
INTEGER :: JI,JJ,JK ! Loop indexes
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZUMEAN_OLD
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZWMEAN_OLD
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZTHMEAN_OLD
!
INTEGER :: IMI !Current model index
!
@@ -103,22 +108,8 @@ IKE=IKU-JPVEXT
!
ZTEMPT = PTHT*(((PPABST)/XP00)**(XRD/XCPD))
!
- XUM_MEAN = PUT + XUM_MEAN
- XVM_MEAN = PVT + XVM_MEAN
- XWM_MEAN = PWT + XWM_MEAN
- XTHM_MEAN = PTHT + XTHM_MEAN
- XTEMPM_MEAN = ZTEMPT + XTEMPM_MEAN
IF (LPASPOL) XSVT_MEAN = PSVT + XSVT_MEAN
IF (CTURB/='NONE') XTKEM_MEAN = PTKET + XTKEM_MEAN
- XPABSM_MEAN = PPABST + XPABSM_MEAN
-!
- XU2_MEAN = PUT**2 + XU2_MEAN
- XV2_MEAN = PVT**2 + XV2_MEAN
- XW2_MEAN = PWT**2 + XW2_MEAN
- XUW_MEAN = PUT*PWT + XUW_MEAN
- XTH2_MEAN = PTHT**2 + XTH2_MEAN
- XTEMP2_MEAN = ZTEMPT**2 + XTEMP2_MEAN
- XPABS2_MEAN = PPABST**2 + XPABS2_MEAN
!
! Wind turbine variables
IMI = GET_CURRENT_MODEL_INDEX()
@@ -137,6 +128,28 @@ IKE=IKU-JPVEXT
!
MEAN_COUNT = MEAN_COUNT + 1
!
+! Save old mean values for covariance computations
+!
+ ZUMEAN_OLD = XUM_MEAN
+ ZWMEAN_OLD = XWM_MEAN
+ ZTHMEAN_OLD = XTHM_MEAN
+! Welford method for variables whom we compute variances
+!
+ CALL UPDATE_WELFORD(MEAN_COUNT,XUM_MEAN,XU2_M2,PUT)
+ CALL UPDATE_WELFORD(MEAN_COUNT,XVM_MEAN,XV2_M2,PVT)
+ CALL UPDATE_WELFORD(MEAN_COUNT,XWM_MEAN,XW2_M2,PWT)
+ CALL UPDATE_WELFORD(MEAN_COUNT,XTHM_MEAN,XTH2_M2,PTHT)
+ CALL UPDATE_WELFORD(MEAN_COUNT,XTEMPM_MEAN,XTEMP2_M2,ZTEMPT)
+ CALL UPDATE_WELFORD(MEAN_COUNT,XPABSM_MEAN,XPABS2_M2,PPABST)
+!
+! Welford method for covariances
+!
+ IF (LCOV_FIELD) THEN
+ XUV_MEAN = XUV_MEAN + (PUT-ZUMEAN_OLD)*(PVT-XVM_MEAN)
+ XUW_MEAN = XUW_MEAN + (PUT-ZUMEAN_OLD)*(PWT-XWM_MEAN)
+ XVW_MEAN = XVW_MEAN + (PWT-ZWMEAN_OLD)*(PVT-XVM_MEAN)
+ XWTH_MEAN = XVW_MEAN + (PWT-ZWMEAN_OLD)*(PTHT-XTHM_MEAN)
+ END IF
!
!-----------------------------------------------------------------------
!
diff --git a/src/MNH/modd_mean_field.f90 b/src/MNH/modd_mean_field.f90
index 8c0e6a0515d9e9cae32b304157e1ba7e8676e30d..71e59e61fc8b88edb54f2785bef8bb05501030a0 100644
--- a/src/MNH/modd_mean_field.f90
+++ b/src/MNH/modd_mean_field.f90
@@ -48,5 +48,6 @@ USE MODD_PARAMETERS
IMPLICIT NONE
LOGICAL :: LMEAN_FIELD
+LOGICAL :: LCOV_FIELD
END MODULE MODD_MEAN_FIELD
diff --git a/src/MNH/modd_mean_fieldn.f90 b/src/MNH/modd_mean_fieldn.f90
index 38572bc0bb2cb8f50ea1bc4cd4352100af779c84..eb346a3e709798b9ee6efe107c7d6b98ca9d1fa0 100644
--- a/src/MNH/modd_mean_fieldn.f90
+++ b/src/MNH/modd_mean_fieldn.f90
@@ -31,6 +31,7 @@
!! -------------
!! Original 01/07/11
!! 10/2016 (C.Lac) Add max values
+!! 11/2022 (E. Jezequel) add covariances
!!
!-------------------------------------------------------------------------------
!
@@ -38,6 +39,7 @@
! ------------
!
USE MODD_PARAMETERS, ONLY: JPMODELMAX
+USE MODD_MEAN_FIELD
IMPLICIT NONE
TYPE MEAN_FIELD_t
@@ -48,10 +50,13 @@ TYPE MEAN_FIELD_t
REAL, DIMENSION(:,:,:), POINTER :: XPABSM_MEAN=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XSVT_MEAN=>NULL()
- REAL, DIMENSION(:,:,:), POINTER :: XU2_MEAN=>NULL(),XV2_MEAN=>NULL(),XW2_MEAN=>NULL(),XUW_MEAN=>NULL()
- REAL, DIMENSION(:,:,:), POINTER :: XTH2_MEAN=>NULL()
- REAL, DIMENSION(:,:,:), POINTER :: XTEMP2_MEAN=>NULL()
- REAL, DIMENSION(:,:,:), POINTER :: XPABS2_MEAN=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XU2_M2=>NULL(),XV2_M2=>NULL(),XW2_M2=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XTH2_M2=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XTEMP2_M2=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XPABS2_M2=>NULL()
+
+ REAL, DIMENSION(:,:,:), POINTER :: XUV_MEAN=>NULL(),XUW_MEAN=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XVW_MEAN=>NULL(),XWTH_MEAN=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XUM_MAX=>NULL(),XVM_MAX=>NULL(),XWM_MAX=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XTHM_MAX=>NULL()
@@ -74,10 +79,13 @@ REAL, DIMENSION(:,:,:), POINTER :: XTKEM_MEAN=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XPABSM_MEAN=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XSVT_MEAN=>NULL()
-REAL, DIMENSION(:,:,:), POINTER :: XU2_MEAN=>NULL(),XV2_MEAN=>NULL(),XW2_MEAN=>NULL(),XUW_MEAN=>NULL()
-REAL, DIMENSION(:,:,:), POINTER :: XTH2_MEAN=>NULL()
-REAL, DIMENSION(:,:,:), POINTER :: XTEMP2_MEAN=>NULL()
-REAL, DIMENSION(:,:,:), POINTER :: XPABS2_MEAN=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XU2_M2=>NULL(),XV2_M2=>NULL(),XW2_M2=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XTH2_M2=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XTEMP2_M2=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XPABS2_M2=>NULL()
+
+REAL, DIMENSION(:,:,:), POINTER :: XUV_MEAN=>NULL(),XUW_MEAN=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XVW_MEAN=>NULL(),XWTH_MEAN=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XUM_MAX=>NULL(),XVM_MAX=>NULL(),XWM_MAX=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XTHM_MAX=>NULL()
@@ -110,14 +118,19 @@ MEAN_FIELD_MODEL(KFROM)%XTEMPM_MAX=>XTEMPM_MAX
MEAN_FIELD_MODEL(KFROM)%XTKEM_MAX=>XTKEM_MAX
MEAN_FIELD_MODEL(KFROM)%XPABSM_MAX=>XPABSM_MAX
-MEAN_FIELD_MODEL(KFROM)%XU2_MEAN=>XU2_MEAN
-MEAN_FIELD_MODEL(KFROM)%XV2_MEAN=>XV2_MEAN
-MEAN_FIELD_MODEL(KFROM)%XW2_MEAN=>XW2_MEAN
-MEAN_FIELD_MODEL(KFROM)%XUW_MEAN=>XUW_MEAN
-MEAN_FIELD_MODEL(KFROM)%XTH2_MEAN=>XTH2_MEAN
-MEAN_FIELD_MODEL(KFROM)%XTEMP2_MEAN=>XTEMP2_MEAN
-MEAN_FIELD_MODEL(KFROM)%XPABS2_MEAN=>XPABS2_MEAN
-
+MEAN_FIELD_MODEL(KFROM)%XU2_M2=>XU2_M2
+MEAN_FIELD_MODEL(KFROM)%XV2_M2=>XV2_M2
+MEAN_FIELD_MODEL(KFROM)%XW2_M2=>XW2_M2
+MEAN_FIELD_MODEL(KFROM)%XTH2_M2=>XTH2_M2
+MEAN_FIELD_MODEL(KFROM)%XTEMP2_M2=>XTEMP2_M2
+MEAN_FIELD_MODEL(KFROM)%XPABS2_M2=>XPABS2_M2
+
+IF (LCOV_FIELD) THEN
+ MEAN_FIELD_MODEL(KFROM)%XUV_MEAN=>XUV_MEAN
+ MEAN_FIELD_MODEL(KFROM)%XUW_MEAN=>XUW_MEAN
+ MEAN_FIELD_MODEL(KFROM)%XVW_MEAN=>XVW_MEAN
+ MEAN_FIELD_MODEL(KFROM)%XWTH_MEAN=>XWTH_MEAN
+END IF
!
! Current model is set to model KTO
XUM_MEAN=>MEAN_FIELD_MODEL(KTO)%XUM_MEAN
@@ -137,13 +150,19 @@ XTEMPM_MAX=>MEAN_FIELD_MODEL(KTO)%XTEMPM_MAX
XTKEM_MAX=>MEAN_FIELD_MODEL(KTO)%XTKEM_MAX
XPABSM_MAX=>MEAN_FIELD_MODEL(KTO)%XPABSM_MAX
-XU2_MEAN=>MEAN_FIELD_MODEL(KTO)%XU2_MEAN
-XV2_MEAN=>MEAN_FIELD_MODEL(KTO)%XV2_MEAN
-XW2_MEAN=>MEAN_FIELD_MODEL(KTO)%XW2_MEAN
-XUW_MEAN=>MEAN_FIELD_MODEL(KTO)%XUW_MEAN
-XTH2_MEAN=>MEAN_FIELD_MODEL(KTO)%XTH2_MEAN
-XTEMP2_MEAN=>MEAN_FIELD_MODEL(KTO)%XTEMP2_MEAN
-XPABS2_MEAN=>MEAN_FIELD_MODEL(KTO)%XPABS2_MEAN
+XU2_M2=>MEAN_FIELD_MODEL(KTO)%XU2_M2
+XV2_M2=>MEAN_FIELD_MODEL(KTO)%XV2_M2
+XW2_M2=>MEAN_FIELD_MODEL(KTO)%XW2_M2
+XTH2_M2=>MEAN_FIELD_MODEL(KTO)%XTH2_M2
+XTEMP2_M2=>MEAN_FIELD_MODEL(KTO)%XTEMP2_M2
+XPABS2_M2=>MEAN_FIELD_MODEL(KTO)%XPABS2_M2
+
+IF (LCOV_FIELD) THEN
+ XUV_MEAN=>MEAN_FIELD_MODEL(KTO)%XUV_MEAN
+ XUW_MEAN=>MEAN_FIELD_MODEL(KTO)%XUW_MEAN
+ XVW_MEAN=>MEAN_FIELD_MODEL(KTO)%XVW_MEAN
+ XWTH_MEAN=>MEAN_FIELD_MODEL(KTO)%XWTH_MEAN
+END IF
MEAN_COUNT=>MEAN_FIELD_MODEL(KTO)%MEAN_COUNT
diff --git a/src/MNH/modn_mean.f90 b/src/MNH/modn_mean.f90
index 99f0143f6c881af2e268642ec9bbe571973788f6..5b39656863e21fafe1030fdeafc3c894a6ae57ca 100644
--- a/src/MNH/modn_mean.f90
+++ b/src/MNH/modn_mean.f90
@@ -38,7 +38,7 @@ USE MODD_MEAN_FIELD
!
IMPLICIT NONE
!
-NAMELIST/NAM_MEAN/LMEAN_FIELD
+NAMELIST/NAM_MEAN/LMEAN_FIELD,LCOV_FIELD
!
!
END MODULE MODN_MEAN
diff --git a/src/MNH/update_welford.f90 b/src/MNH/update_welford.f90
new file mode 100644
index 0000000000000000000000000000000000000000..8aa8e869b090caebb1d305676188606df19c20a5
--- /dev/null
+++ b/src/MNH/update_welford.f90
@@ -0,0 +1,63 @@
+!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.
+!-----------------------------------------------------------------
+! #######################
+ MODULE MODI_UPDATE_WELFORD
+! #######################
+!
+INTERFACE
+!
+SUBROUTINE UPDATE_WELFORD(KCOUNT,PMEAN,PM2,PVALUE)
+ INTEGER, INTENT(IN) :: KCOUNT ! time count
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PMEAN ! aggregate for mean value
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PM2 ! aggregate for var value
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PVALUE ! variables
+END SUBROUTINE UPDATE_WELFORD
+!
+END INTERFACE
+!
+END MODULE MODI_UPDATE_WELFORD
+!#########################################################
+SUBROUTINE UPDATE_WELFORD(KCOUNT,PMEAN,PM2,PVALUE)
+!!**** *UPDATE_WELFORD* - Welford Algorithm
+!!
+!! PURPOSE
+!! -------
+! The purpose of this module is to compute variance on a large number of sample
+! and avoid catastrophic cancellation, leading to negative variance
+!
+!!
+!!** IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!! E. Jezequel
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original November, 2022
+! ------------------------------------------------------------------------------
+!
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KCOUNT ! time count
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PMEAN ! aggregate for mean value
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PM2 ! aggregate for var value
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PVALUE ! current value
+! local variables
+REAL, DIMENSION(SIZE(PVALUE,1),SIZE(PVALUE,2),SIZE(PVALUE,3)) :: ZDELTA, ZDELTA2
+!
+ZDELTA = PVALUE - PMEAN
+PMEAN = PMEAN + ZDELTA / KCOUNT
+ZDELTA2= PVALUE - PMEAN
+PM2 = PM2 + ZDELTA * ZDELTA2
+!
+END SUBROUTINE UPDATE_WELFORD
+!
diff --git a/src/MNH/write_lfin.f90 b/src/MNH/write_lfin.f90
index 4e11a572ddfe69d8ea56805d977068321fd4efbb..dd96607f1733505d37f5e5a40a110a1f5729cec1 100644
--- a/src/MNH/write_lfin.f90
+++ b/src/MNH/write_lfin.f90
@@ -181,6 +181,7 @@ END MODULE MODI_WRITE_LFIFM_n
! P. Wautelet 11/03/2021: bugfix: correct name for NSV_LIMA_IMM_NUCL
! J.L. Redelsperger 03/2021: add OCEAN and auto-coupled O-A LES cases
! A. Costes 12/2021: add Blaze fire model
+! E. Jezequel 11/2022 : add covariances from MEAN fields
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -195,6 +196,7 @@ USE MODD_GRID_n
USE MODD_TIME
USE MODD_TIME_n
USE MODD_FIELD_n
+USE MODD_MEAN_FIELD
USE MODD_MEAN_FIELD_n
USE MODD_DUMMY_GR_FIELD_n
USE MODD_LSFIELD_n
@@ -722,14 +724,14 @@ IF (MEAN_COUNT /= 0) THEN
TZFIELD%CLONGNAME = 'UMME'
TZFIELD%CUNITS = 'm s-1'
TZFIELD%CCOMMENT = 'X_Y_Z_U component of mean wind'
- ZWORK3D = XUM_MEAN/MEAN_COUNT
+ ZWORK3D = XUM_MEAN
CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'U2ME'
TZFIELD%CLONGNAME = 'U2ME'
TZFIELD%CUNITS = 'm2 s-2'
TZFIELD%CCOMMENT = 'X_Y_Z_U component of mean wind variance'
- ZWORK3D = XU2_MEAN/MEAN_COUNT-XUM_MEAN**2/MEAN_COUNT**2
+ ZWORK3D = XU2_M2/MEAN_COUNT
CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'UMMA'
@@ -738,27 +740,51 @@ IF (MEAN_COUNT /= 0) THEN
TZFIELD%CCOMMENT = 'X_Y_Z_U component of max wind'
CALL IO_Field_write(TPFILE,TZFIELD,XUM_MAX)
!
- TZFIELD%CMNHNAME = 'UWME'
- TZFIELD%CLONGNAME = 'UWME'
- TZFIELD%CUNITS = 'm2 s-2'
- TZFIELD%CCOMMENT = 'X_Y_Z_UW component of mean wind variance'
- ZWORK3D = XUW_MEAN/MEAN_COUNT-(XUM_MEAN*XWM_MEAN)/MEAN_COUNT**2
- CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
- !
+ IF (LCOV_FIELD) THEN
+ !
+ TZFIELD%CMNHNAME = 'UVME'
+ TZFIELD%CLONGNAME = 'UVME'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CCOMMENT = 'X_Y_Z_UV component of mean wind variance'
+ ZWORK3D = XUV_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+ !
+ TZFIELD%CMNHNAME = 'UWME'
+ TZFIELD%CLONGNAME = 'UWME'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CCOMMENT = 'X_Y_Z_UW component of mean wind variance'
+ ZWORK3D = XUW_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+ !
+ TZFIELD%CMNHNAME = 'VWME'
+ TZFIELD%CLONGNAME = 'VWME'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CCOMMENT = 'X_Y_Z_VW component of mean wind variance'
+ ZWORK3D = XVW_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+ !
+ TZFIELD%CMNHNAME = 'WTHME'
+ TZFIELD%CLONGNAME = 'WTHME'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CCOMMENT = 'X_Y_Z_WTH component of mean wind variance'
+ ZWORK3D = XWTH_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+ !
+ END IF
TZFIELD%NGRID = 3
!
TZFIELD%CMNHNAME = 'VMME'
TZFIELD%CLONGNAME = 'VMME'
TZFIELD%CUNITS = 'm s-1'
TZFIELD%CCOMMENT = 'X_Y_Z_V component of mean wind'
- ZWORK3D = XVM_MEAN/MEAN_COUNT
+ ZWORK3D = XVM_MEAN
CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'V2ME'
TZFIELD%CLONGNAME = 'V2ME'
TZFIELD%CUNITS = 'm2 s-2'
TZFIELD%CCOMMENT = 'X_Y_Z_V component of mean wind variance'
- ZWORK3D = XV2_MEAN/MEAN_COUNT-XVM_MEAN**2/MEAN_COUNT**2
+ ZWORK3D = XV2_M2/MEAN_COUNT
CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'VMMA'
@@ -773,14 +799,14 @@ IF (MEAN_COUNT /= 0) THEN
TZFIELD%CLONGNAME = 'WMME'
TZFIELD%CUNITS = 'm s-1'
TZFIELD%CCOMMENT = 'X_Y_Z_vertical mean wind'
- ZWORK3D = XWM_MEAN/MEAN_COUNT
+ ZWORK3D = XWM_MEAN
CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'W2ME'
TZFIELD%CLONGNAME = 'W2ME'
TZFIELD%CUNITS = 'm2 s-2'
TZFIELD%CCOMMENT = 'X_Y_Z_vertical mean wind variance'
- ZWORK3D = XW2_MEAN/MEAN_COUNT-XWM_MEAN**2/MEAN_COUNT**2
+ ZWORK3D = XW2_M2/MEAN_COUNT
CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'WMMA'
@@ -802,14 +828,14 @@ IF (MEAN_COUNT /= 0) THEN
TZFIELD%CLONGNAME = 'THMME'
TZFIELD%CUNITS = 'K'
TZFIELD%CCOMMENT = 'X_Y_Z_mean potential temperature'
- ZWORK3D = XTHM_MEAN/MEAN_COUNT
+ ZWORK3D = XTHM_MEAN
CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'TH2ME'
TZFIELD%CLONGNAME = 'TH2ME'
TZFIELD%CUNITS = 'K2'
TZFIELD%CCOMMENT = 'X_Y_Z_mean potential temperature variance'
- ZWORK3D = XTH2_MEAN/MEAN_COUNT-XTHM_MEAN**2/MEAN_COUNT**2
+ ZWORK3D = XTH2_M2/MEAN_COUNT
CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'THMMA'
@@ -822,14 +848,14 @@ IF (MEAN_COUNT /= 0) THEN
TZFIELD%CLONGNAME = 'TEMPMME'
TZFIELD%CUNITS = 'K'
TZFIELD%CCOMMENT = 'X_Y_Z_mean temperature'
- ZWORK3D = XTEMPM_MEAN/MEAN_COUNT
+ ZWORK3D = XTEMPM_MEAN
CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'TEMP2ME'
TZFIELD%CLONGNAME = 'TEMP2ME'
TZFIELD%CUNITS = 'K2'
TZFIELD%CCOMMENT = 'X_Y_Z_mean temperature variance'
- ZWORK3D = XTEMP2_MEAN/MEAN_COUNT-XTEMPM_MEAN**2/MEAN_COUNT**2
+ ZWORK3D = XTEMP2_M2/MEAN_COUNT
CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'TEMPMMA'
@@ -842,14 +868,14 @@ IF (MEAN_COUNT /= 0) THEN
TZFIELD%CLONGNAME = 'PABSMME'
TZFIELD%CUNITS = 'Pa'
TZFIELD%CCOMMENT = 'X_Y_Z_mean ABSolute Pressure'
- ZWORK3D = XPABSM_MEAN/MEAN_COUNT
+ ZWORK3D = XPABSM_MEAN
CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'PABS2ME'
TZFIELD%CLONGNAME = 'PABS2ME'
TZFIELD%CUNITS = 'Pa2'
TZFIELD%CCOMMENT = 'X_Y_Z_mean ABSolute Pressure variance'
- ZWORK3D = XPABS2_MEAN/MEAN_COUNT-XPABSM_MEAN**2/MEAN_COUNT**2
+ ZWORK3D = XPABS2_M2/MEAN_COUNT
CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'PABSMMA'