From f16c676ffea6088a93d2e8269b3c35f434f1fbb5 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Tue, 22 Feb 2022 11:51:54 +0100
Subject: [PATCH] Quentin 22/02/2022: add turb files from common and aux from
MNH-55 before modifs
---
src/mesonh/aux/gradient_m.f90 | 744 +++++++++++++++
src/mesonh/aux/gradient_u.f90 | 323 +++++++
src/mesonh/aux/gradient_v.f90 | 321 +++++++
src/mesonh/aux/gradient_w.f90 | 297 ++++++
src/mesonh/turb/mode_turb_ver_dyn_flux.f90 | 778 ++++++++++++++++
src/mesonh/turb/mode_turb_ver_sv_corr.f90 | 178 ++++
src/mesonh/turb/mode_turb_ver_sv_flux.f90 | 446 +++++++++
src/mesonh/turb/mode_turb_ver_thermo_corr.f90 | 866 ++++++++++++++++++
src/mesonh/turb/mode_turb_ver_thermo_flux.f90 | 833 +++++++++++++++++
9 files changed, 4786 insertions(+)
create mode 100644 src/mesonh/aux/gradient_m.f90
create mode 100644 src/mesonh/aux/gradient_u.f90
create mode 100644 src/mesonh/aux/gradient_v.f90
create mode 100644 src/mesonh/aux/gradient_w.f90
create mode 100644 src/mesonh/turb/mode_turb_ver_dyn_flux.f90
create mode 100644 src/mesonh/turb/mode_turb_ver_sv_corr.f90
create mode 100644 src/mesonh/turb/mode_turb_ver_sv_flux.f90
create mode 100644 src/mesonh/turb/mode_turb_ver_thermo_corr.f90
create mode 100644 src/mesonh/turb/mode_turb_ver_thermo_flux.f90
diff --git a/src/mesonh/aux/gradient_m.f90 b/src/mesonh/aux/gradient_m.f90
new file mode 100644
index 000000000..b5ec025aa
--- /dev/null
+++ b/src/mesonh/aux/gradient_m.f90
@@ -0,0 +1,744 @@
+!MNH_LIC Copyright 1994-2020 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_GRADIENT_M
+! ######################
+!
+INTERFACE
+!
+!
+FUNCTION GX_M_M(PA,PDXX,PDZZ,PDZX) RESULT(PGX_M_M)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_M_M ! result mass point
+!
+END FUNCTION GX_M_M
+!
+!
+FUNCTION GY_M_M(PA,PDYY,PDZZ,PDZY) RESULT(PGY_M_M)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_M_M ! result mass point
+!
+END FUNCTION GY_M_M
+!
+!
+FUNCTION GZ_M_M(PA,PDZZ) RESULT(PGZ_M_M)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGZ_M_M ! result mass point
+!
+END FUNCTION GZ_M_M
+!
+ FUNCTION GX_M_U(KKA,KKU,KL,PY,PDXX,PDZZ,PDZX) RESULT(PGX_M_U)
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN) :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGX_M_U ! result at flux
+ ! side
+END FUNCTION GX_M_U
+!
+!
+ FUNCTION GY_M_V(KKA,KKU,KL,PY,PDYY,PDZZ,PDZY) RESULT(PGY_M_V)
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN) :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY !d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY !d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ !d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGY_M_V ! result at flux
+ ! side
+END FUNCTION GY_M_V
+!
+ FUNCTION GZ_M_W(KKA,KKU,KL,PY,PDZZ) RESULT(PGZ_M_W)
+!
+IMPLICIT NONE
+!
+ ! Metric coefficient:
+INTEGER, INTENT(IN) :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN) :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ !d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGZ_M_W ! result at flux
+ ! side
+!
+END FUNCTION GZ_M_W
+!
+END INTERFACE
+!
+END MODULE MODI_GRADIENT_M
+!
+!
+!
+! #######################################################
+ FUNCTION GX_M_M(PA,PDXX,PDZZ,PDZX) RESULT(PGX_M_M)
+! #######################################################
+!
+!!**** *GX_M_M* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian X
+!! direction for a variable placed at the
+!! mass point and the result is placed at
+!! the mass point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the X cartesian direction for a field PA placed at the
+! mass point. The result is placed at the mass point.
+!
+!
+! ( ______________z )
+! ( (___x ) )
+! 1 ( _x (d*zx dzm(PA) ) )
+! PGX_M_M = ---- (dxf(PA) - (------------)) )
+! ___x ( ( ) )
+! d*xx ( ( d*zz ) )
+!
+!
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MXM,MXF,MZF : Shuman functions (mean operators)
+!! DXF,DZF : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! MODD_CONF : LFLAT
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 18/07/94
+!! 19/07/00 add the LFLAT switch (J. Stein)
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+USE MODD_CONF, ONLY:LFLAT
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_M_M ! result mass point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GX_M_M
+! --------------------
+!
+IF (.NOT. LFLAT) THEN
+ PGX_M_M(:,:,:)= (DXF(MXM(PA(:,:,:))) - &
+ MZF(MXF(PDZX)*DZM(PA(:,:,:)) &
+ /PDZZ(:,:,:)) ) /MXF(PDXX(:,:,:))
+ELSE
+ PGX_M_M(:,:,:)=DXF(MXM(PA(:,:,:))) / MXF(PDXX(:,:,:))
+END IF
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GX_M_M
+!
+!
+! #######################################################
+ FUNCTION GY_M_M(PA,PDYY,PDZZ,PDZY) RESULT(PGY_M_M)
+! #######################################################
+!
+!!**** *GY_M_M* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Y
+!! direction for a variable placed at the
+!! mass point and the result is placed at
+!! the mass point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the Y cartesian direction for a field PA placed at the
+! mass point. The result is placed at the mass point.
+!
+!
+! ( ______________z )
+! ( (___y ) )
+! 1 ( _y (d*zy dzm(PA) ) )
+! PGY_M_M = ---- (dyf(PA) - (------------)) )
+! ___y ( ( ) )
+! d*yy ( ( d*zz ) )
+!
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MYM,MYF,MZF : Shuman functions (mean operators)
+!! DYF,DZF : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! MODD_CONF : LFLAT
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 18/07/94
+!! 19/07/00 add the LFLAT switch (J. Stein)
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODD_CONF, ONLY:LFLAT
+USE MODI_SHUMAN
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_M_M ! result mass point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GY_M_M
+! --------------------
+!
+!
+IF (.NOT. LFLAT) THEN
+ PGY_M_M(:,:,:)= (DYF(MYM(PA))-MZF(MYF(PDZY)*DZM(PA)&
+ /PDZZ) ) /MYF(PDYY)
+ELSE
+ PGY_M_M(:,:,:)= DYF(MYM(PA))/MYF(PDYY)
+ENDIF
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GY_M_M
+
+!
+!
+!
+! #############################################
+ FUNCTION GZ_M_M(PA,PDZZ) RESULT(PGZ_M_M)
+! #############################################
+!
+!!**** *GZ_M_M* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Z
+!! direction for a variable placed at the
+!! mass point and the result is placed at
+!! the mass point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the Z cartesian direction for a field PA placed at the
+! mass point. The result is placed at the mass point.
+!
+! _________z
+! (dzm(PA))
+! PGZ_M_M = (------ )
+! ( d*zz )
+!
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MZF : Shuman functions (mean operators)
+!! DZM : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 18/07/94
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the mass point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGZ_M_M ! result mass point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GZ_M_M
+! --------------------
+!
+PGZ_M_M(:,:,:)= MZF( DZM(PA(:,:,:))/PDZZ(:,:,:) )
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GZ_M_M
+!
+!
+! ##################################################
+ FUNCTION GX_M_U(KKA,KKU,KL,PY,PDXX,PDZZ,PDZX) RESULT(PGX_M_U)
+! ##################################################
+!
+!!**** *GX_M_U * - Compute the gradient along x for a variable localized at
+!! a mass point
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute a gradient along x
+! direction for a field PY localized at a mass point. The result PGX_M_U
+! is localized at a x-flux point (u point).
+!
+! ( ____________z )
+! ( ________x )
+! 1 ( dzm(PY) )
+! PGX_M_U = ---- (dxm(PY) - d*zx -------- )
+! d*xx ( d*zz )
+!
+!
+!
+!!** METHOD
+!! ------
+!! We employ the Shuman operators to compute the derivatives and the
+!! averages. The metric coefficients PDXX,PDZX,PDZZ are dummy arguments.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! FUNCTION DXM: compute a finite difference along the x direction for
+!! a variable at a mass localization
+!! FUNCTION DZM: compute a finite difference along the y direction for
+!! a variable at a mass localization
+!! FUNCTION MXM: compute an average in the x direction for a variable
+!! at a mass localization
+!! FUNCTION MZF: compute an average in the z direction for a variable
+!! at a flux side
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! MODD_CONF : LFLAT
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (function GX_M_U)
+!!
+!!
+!! AUTHOR
+!! ------
+!! P. Hereil and J. Stein * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 05/07/94
+!! Modification 16/03/95 change the order of the arguments
+!! 19/07/00 add the LFLAT switch + inlining(J. Stein)
+!! 20/08/00 optimization (J. Escobar)
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODI_SHUMAN
+USE MODD_CONF, ONLY:LFLAT
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and result
+! ------------------------------------
+!
+INTEGER, INTENT(IN) :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN) :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGX_M_U ! result at flux
+ ! side
+INTEGER IIU,IKU,JI,JK
+!
+INTEGER :: JJK,IJU
+INTEGER :: JIJK,JIJKOR,JIJKEND
+INTEGER :: JI_1JK, JIJK_1, JI_1JK_1, JIJKP1, JI_1JKP1
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE THE GRADIENT ALONG X
+! -----------------------------
+!
+IIU=SIZE(PY,1)
+IJU=SIZE(PY,2)
+IKU=SIZE(PY,3)
+IF (.NOT. LFLAT) THEN
+! PGX_M_U = ( DXM(PY) - MZF ( MXM( DZM(PY) /PDZZ ) * PDZX ) )/PDXX
+!! DO JK=1+JPVEXT_TURB,IKU-JPVEXT_TURB
+!! DO JI=1+JPHEXT,IIU
+!! PGX_M_U(JI,:,JK)= &
+!! ( PY(JI,:,JK)-PY(JI-1,:,JK) &
+!! -( (PY(JI,:,JK)-PY(JI,:,JK-1)) / PDZZ(JI,:,JK) &
+!! +(PY(JI-1,:,JK)-PY(JI-1,:,JK-1)) / PDZZ(JI-1,:,JK) &
+!! ) * PDZX(JI,:,JK)* 0.25 &
+!! -( (PY(JI,:,JK+1)-PY(JI,:,JK)) / PDZZ(JI,:,JK+1) &
+!! +(PY(JI-1,:,JK+1)-PY(JI-1,:,JK)) / PDZZ(JI-1,:,JK+1) &
+!! ) * PDZX(JI,:,JK+1)* 0.25 &
+!! ) / PDXX(JI,:,JK)
+!! END DO
+!! END DO
+ JIJKOR = 1 + JPHEXT + IIU*IJU*(JPVEXT_TURB+1 - 1)
+ JIJKEND = IIU*IJU*(IKU-JPVEXT_TURB)
+!CDIR NODEP
+!OCL NOVREC
+ DO JIJK=JIJKOR , JIJKEND
+! indexation
+ JI_1JK = JIJK - 1
+ JIJK_1 = JIJK - IIU*IJU*KL
+ JI_1JK_1 = JIJK - 1 - IIU*IJU*KL
+ JIJKP1 = JIJK + IIU*IJU*KL
+ JI_1JKP1 = JIJK - 1 + IIU*IJU*KL
+!
+ PGX_M_U(JIJK,1,1)= &
+ ( PY(JIJK,1,1)-PY(JI_1JK,1,1) &
+ -( (PY(JIJK,1,1)-PY(JIJK_1,1,1)) / PDZZ(JIJK,1,1) &
+ +(PY(JI_1JK,1,1)-PY(JI_1JK_1,1,1)) / PDZZ(JI_1JK,1,1) &
+ ) * PDZX(JIJK,1,1)* 0.25 &
+ -( (PY(JIJKP1,1,1)-PY(JIJK,1,1)) / PDZZ(JIJKP1,1,1) &
+ +(PY(JI_1JKP1,1,1)-PY(JI_1JK,1,1)) / PDZZ(JI_1JKP1,1,1) &
+ ) * PDZX(JIJKP1,1,1)* 0.25 &
+ ) / PDXX(JIJK,1,1)
+ END DO
+
+!
+ DO JI=1+JPHEXT,IIU
+ PGX_M_U(JI,:,KKU)= ( PY(JI,:,KKU)-PY(JI-1,:,KKU) ) / PDXX(JI,:,KKU)
+ PGX_M_U(JI,:,KKA)= PGX_M_U(JI,:,KKU) ! -999.
+ END DO
+ELSE
+! PGX_M_U = DXM(PY) / PDXX
+ PGX_M_U(1+1:IIU,:,:) = ( PY(1+1:IIU,:,:)-PY(1:IIU-1,:,:) ) & ! +JPHEXT
+ / PDXX(1+1:IIU,:,:)
+ENDIF
+DO JI=1,JPHEXT
+ PGX_M_U(JI,:,:)=PGX_M_U(IIU-2*JPHEXT+JI,:,:) ! for reprod JPHEXT <> 1
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION GX_M_U
+!
+!
+! ##################################################
+ FUNCTION GY_M_V(KKA,KKU,KL,PY,PDYY,PDZZ,PDZY) RESULT(PGY_M_V)
+! ##################################################
+!
+!!**** *GY_M_V * - Compute the gradient along y for a variable localized at
+!! a mass point
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute a gradient along y
+! direction for a field PY localized at a mass point. The result PGY_M_V
+! is localized at a y-flux point (v point).
+!
+! ( ____________z )
+! ( ________y )
+! 1 ( dzm(PY) )
+! PGY_M_V = ---- (dym(PY) - d*zy -------- )
+! d*yy ( d*zz )
+!
+!
+!
+!
+!!** METHOD
+!! ------
+!! We employ the Shuman operators to compute the derivatives and the
+!! averages. The metric coefficients PDYY,PDZY,PDZZ are dummy arguments.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! FUNCTION DYM: compute a finite difference along the y direction for
+!! a variable at a mass localization
+!! FUNCTION DZM: compute a finite difference along the y direction for
+!! a variable at a mass localization
+!! FUNCTION MYM: compute an average in the x direction for a variable
+!! at a mass localization
+!! FUNCTION MZF: compute an average in the z direction for a variable
+!! at a flux side
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! MODD_CONF : LFLAT
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (function GY_M_V)
+!!
+!!
+!! AUTHOR
+!! ------
+!! P. Hereil and J. Stein * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 05/07/94
+!! Modification 16/03/95 change the order of the arguments
+!! 19/07/00 add the LFLAT switch + inlining(J. Stein)
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODI_SHUMAN
+USE MODD_CONF, ONLY:LFLAT
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+! -------------------------------------
+!
+INTEGER, INTENT(IN) :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN) :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY !d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY !d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ !d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGY_M_V ! result at flux
+ ! side
+INTEGER IJU,IKU,JJ,JK
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE THE GRADIENT ALONG Y
+! ----------------------------
+!
+IJU=SIZE(PY,2)
+IKU=SIZE(PY,3)
+IF (.NOT. LFLAT) THEN
+! PGY_M_V = ( DYM(PY) - MZF ( MYM( DZM(PY) /PDZZ ) * PDZY ) )/PDYY
+ DO JK=1+JPVEXT_TURB,IKU-JPVEXT_TURB
+ DO JJ=1+JPHEXT,IJU
+ PGY_M_V(:,JJ,JK)= &
+ ( PY(:,JJ,JK)-PY(:,JJ-1,JK) &
+ -( (PY(:,JJ,JK)-PY(:,JJ,JK-KL)) / PDZZ(:,JJ,JK) &
+ +(PY(:,JJ-1,JK)-PY(:,JJ-KL,JK-KL)) / PDZZ(:,JJ-1,JK) &
+ ) * PDZY(:,JJ,JK)* 0.25 &
+ -( (PY(:,JJ,JK+KL)-PY(:,JJ,JK)) / PDZZ(:,JJ,JK+KL) &
+ +(PY(:,JJ-1,JK+KL)-PY(:,JJ-1,JK)) / PDZZ(:,JJ-1,JK+KL) &
+ ) * PDZY(:,JJ,JK+KL)* 0.25 &
+ ) / PDYY(:,JJ,JK)
+ END DO
+ END DO
+!
+ DO JJ=1+JPHEXT,IJU
+ PGY_M_V(:,JJ,KKU)= ( PY(:,JJ,KKU)-PY(:,JJ-1,KKU) ) / PDYY(:,JJ,KKU)
+ PGY_M_V(:,JJ,KKA)= PGY_M_V(:,JJ,KKU) ! -999.
+ END DO
+ELSE
+! PGY_M_V = DYM(PY)/PDYY
+ PGY_M_V(:,1+1:IJU,:) = ( PY(:,1+1:IJU,:)-PY(:,1:IJU-1,:) ) & ! +JPHEXT
+ / PDYY(:,1+1:IJU,:)
+ENDIF
+DO JJ=1,JPHEXT
+ PGY_M_V(:,JJ,:)=PGY_M_V(:,IJU-2*JPHEXT+JJ,:)
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION GY_M_V
+!
+!
+! #########################################
+ FUNCTION GZ_M_W(KKA,KKU,KL,PY,PDZZ) RESULT(PGZ_M_W)
+! #########################################
+!
+!!**** *GZ_M_W * - Compute the gradient along z direction for a
+!! variable localized at a mass point
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute a gradient along x,y,z
+! directions for a field PY localized at a mass point. The result PGZ_M_W
+! is localized at a z-flux point (w point)
+!
+!
+! dzm(PY)
+! PGZ_M_W = -------
+! d*zz
+!
+!!** METHOD
+!! ------
+!! We employ the Shuman operators to compute the derivatives and the
+!! averages. The metric coefficients PDZZ are dummy arguments.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! FUNCTION DZM : compute a finite difference along the z
+!! direction for a variable at a mass localization
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODI_SHUMAN : interface for the Shuman functions
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (function GZ_M_W)
+!!
+!!
+!! AUTHOR
+!! ------
+!! P. Hereil and J. Stein * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 05/07/94
+!! Modification 16/03/95 change the order of the arguments
+!! 19/07/00 inlining(J. Stein)
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODI_SHUMAN
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments and results
+! -------------------------------------
+!
+INTEGER, INTENT(IN) :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN) :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
+
+ ! Metric coefficient:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ !d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! variable at mass
+ ! localization
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGZ_M_W ! result at flux
+ ! side
+!
+INTEGER :: IKT,IKTB,IKTE
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE THE GRADIENT ALONG Z
+! -----------------------------
+!
+IKT=SIZE(PY,3)
+IKTB=1+JPVEXT_TURB
+IKTE=IKT-JPVEXT_TURB
+
+PGZ_M_W(:,:,IKTB:IKTE) = (PY(:,:,IKTB:IKTE)-PY(:,:,IKTB-KL:IKTE-KL)) &
+ / PDZZ(:,:,IKTB:IKTE)
+PGZ_M_W(:,:,KKU)= (PY(:,:,KKU)-PY(:,:,KKU-KL)) &
+ / PDZZ(:,:,KKU)
+PGZ_M_W(:,:,KKA)= PGZ_M_W(:,:,KKU) ! -999.
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION GZ_M_W
+
diff --git a/src/mesonh/aux/gradient_u.f90 b/src/mesonh/aux/gradient_u.f90
new file mode 100644
index 000000000..3d32ffa80
--- /dev/null
+++ b/src/mesonh/aux/gradient_u.f90
@@ -0,0 +1,323 @@
+!MNH_LIC Copyright 1994-2020 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_GRADIENT_U
+! ######################
+!
+INTERFACE
+!
+!
+FUNCTION GX_U_M(PA,PDXX,PDZZ,PDZX) RESULT(PGX_U_M)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the U point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_U_M ! result mass point
+!
+END FUNCTION GX_U_M
+!
+!
+FUNCTION GY_U_UV(PA,PDYY,PDZZ,PDZY) RESULT(PGY_U_UV)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the U point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_U_UV ! result UV point
+!
+END FUNCTION GY_U_UV
+!
+!
+FUNCTION GZ_U_UW(PA,PDZZ) RESULT(PGZ_U_UW)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the U point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGZ_U_UW ! result UW point
+!
+END FUNCTION GZ_U_UW
+!
+END INTERFACE
+!
+END MODULE MODI_GRADIENT_U
+!
+!
+!
+!
+! #######################################################
+ FUNCTION GX_U_M(PA,PDXX,PDZZ,PDZX) RESULT(PGX_U_M)
+! #######################################################
+!
+!!**** *GX_U_M* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian X
+!! direction for a variable placed at the
+!! U point and the result is placed at
+!! the mass point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the X cartesian direction for a field PA placed at the
+! U point. The result is placed at the mass point.
+!
+!
+! ( ______________z )
+! ( (___________x ) )
+! 1 ( (d*zx dzm(PA) ) )
+! PGX_U_M = ---- (dxf(PA) - (------------)) )
+! ___x ( ( ) )
+! d*xx ( ( d*zz ) )
+!
+!
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MXF,MZF : Shuman functions (mean operators)
+!! DXF,DZF : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 19/07/94
+!! 18/10/00 (V.Masson) add LFLAT switch
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+USE MODD_CONF
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the U point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_U_M ! result mass point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GX_U_M
+! --------------------
+!
+IF (.NOT. LFLAT) THEN
+ PGX_U_M(:,:,:)= ( DXF(PA) - &
+ MZF(MXF(PDZX*DZM(PA)) / PDZZ ) &
+ ) / MXF(PDXX)
+ELSE
+ PGX_U_M(:,:,:)= DXF(PA) / MXF(PDXX)
+END IF
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GX_U_M
+!
+!
+! #########################################################
+ FUNCTION GY_U_UV(PA,PDYY,PDZZ,PDZY) RESULT(PGY_U_UV)
+! #########################################################
+!
+!!**** *GY_U_UV* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Y
+!! direction for a variable placed at the
+!! U point and the result is placed at
+!! the UV vorticity point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the Y cartesian direction for a field PA placed at the
+! U point. The result is placed at the UV vorticity point.
+!
+!
+!
+! ( _________________z )
+! ( (___x _________y ) )
+! 1 ( (d*zy (dzm(PA))) ) )
+! PGY_U_UV= ---- (dym(PA) - ( (------ ) ) )
+! ___x ( ( ( ___x ) ) )
+! d*yy ( ( ( d*zz ) ) )
+!
+!
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MXM,MYM,MZF : Shuman functions (mean operators)
+!! DYM,DZM : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 20/07/94
+!! 18/10/00 (V.Masson) add LFLAT switch
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+USE MODD_CONF
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the U point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_U_UV ! result UV point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GY_U_UV
+! ---------------------
+!
+IF (.NOT. LFLAT) THEN
+ PGY_U_UV(:,:,:)= (DYM(PA)- MZF( MYM( DZM(PA)/&
+ MXM(PDZZ) ) *MXM(PDZY) ) ) / MXM(PDYY)
+ELSE
+ PGY_U_UV(:,:,:)= DYM(PA) / MXM(PDYY)
+END IF
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GY_U_UV
+!
+!
+! #######################################################
+ FUNCTION GZ_U_UW(PA,PDZZ) RESULT(PGZ_U_UW)
+! #######################################################
+!
+!!**** *GZ_U_UW - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Z
+!! direction for a variable placed at the
+!! U point and the result is placed at
+!! the UW vorticity point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the Z cartesian direction for a field PA placed at the
+! U point. The result is placed at the UW vorticity point.
+!
+! dzm(PA)
+! PGZ_U_UW = ------
+! ____x
+! d*zz
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MXM : Shuman functions (mean operators)
+!! DZM : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 20/07/94
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the U point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGZ_U_UW ! result UW point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GZ_U_UW
+! ---------------------
+!
+PGZ_U_UW(:,:,:)= DZM(PA) / MXM(PDZZ)
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GZ_U_UW
diff --git a/src/mesonh/aux/gradient_v.f90 b/src/mesonh/aux/gradient_v.f90
new file mode 100644
index 000000000..12c1be749
--- /dev/null
+++ b/src/mesonh/aux/gradient_v.f90
@@ -0,0 +1,321 @@
+!MNH_LIC Copyright 1994-2020 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_GRADIENT_V
+! ######################
+!
+INTERFACE
+!
+!
+FUNCTION GY_V_M(PA,PDYY,PDZZ,PDZY) RESULT(PGY_V_M)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the V point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_V_M ! result mass point
+!
+END FUNCTION GY_V_M
+!
+FUNCTION GX_V_UV(PA,PDXX,PDZZ,PDZX) RESULT(PGX_V_UV)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the V point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_V_UV ! result UV point
+!
+END FUNCTION GX_V_UV
+!
+!
+FUNCTION GZ_V_VW(PA,PDZZ) RESULT(PGZ_V_VW)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the V point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGZ_V_VW ! result VW point
+!
+END FUNCTION GZ_V_VW
+!
+!
+END INTERFACE
+!
+END MODULE MODI_GRADIENT_V
+!
+!
+!
+! #######################################################
+ FUNCTION GY_V_M(PA,PDYY,PDZZ,PDZY) RESULT(PGY_V_M)
+! #######################################################
+!
+!!**** *GY_V_M* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Y
+!! direction for a variable placed at the
+!! V point and the result is placed at
+!! the mass point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the Y cartesian direction for a field PA placed at the
+! V point. The result is placed at the mass point.
+!
+!
+! ( ______________z )
+! ( (___________y ) )
+! 1 ( (d*zy dzm(PA) ) )
+! PGY_V_M = ---- (dyf(PA) - (------------)) )
+! ___y ( ( ) )
+! d*yy ( ( d*zz ) )
+!
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MYF,MZF : Shuman functions (mean operators)
+!! DYF,DZF : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 19/07/94
+!! 18/10/00 (V.Masson) add LFLAT switch
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+USE MODD_CONF
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the V point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_V_M ! result mass point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GY_V_M
+! --------------------
+!
+IF (.NOT. LFLAT) THEN
+ PGY_V_M(:,:,:)= (DYF(PA) - &
+ MZF( MYF(PDZY*DZM(PA))/PDZZ ) &
+ ) / MYF(PDYY)
+ELSE
+ PGY_V_M(:,:,:)= DYF(PA) / MYF(PDYY)
+END IF
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GY_V_M
+!
+!
+! #########################################################
+ FUNCTION GX_V_UV(PA,PDXX,PDZZ,PDZX) RESULT(PGX_V_UV)
+! #########################################################
+!
+!!**** *GX_V_UV* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian X
+!! direction for a variable placed at the
+!! V point and the result is placed at
+!! the UV vorticity point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the X cartesian direction for a field PA placed at the
+! V point. The result is placed at the UV vorticity point.
+!
+!
+! ( _________________z )
+! ( (___y _________x ) )
+! 1 ( (d*zx (dzm(PA))) ) )
+! PGX_V_UV= ---- (dxm(PA) - ( (------ ) ) )
+! ___y ( ( ( ___y ) ) )
+! d*xx ( ( ( d*zz ) ) )
+!
+!
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MXM,MZF,MYM : Shuman functions (mean operators)
+!! DXM,DZM : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 20/07/94
+!! 18/10/00 (V.Masson) add LFLAT switch
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+USE MODD_CONF
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the V point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_V_UV ! result UV point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GX_V_UV
+! ---------------------
+!
+IF (.NOT. LFLAT) THEN
+ PGX_V_UV(:,:,:)= ( DXM(PA)- MZF( MXM( DZM(PA)/&
+ MYM(PDZZ) ) *MYM(PDZX) ) ) / MYM(PDXX)
+ELSE
+ PGX_V_UV(:,:,:)= DXM(PA) / MYM(PDXX)
+END IF
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GX_V_UV
+!
+!
+! #######################################################
+ FUNCTION GZ_V_VW(PA,PDZZ) RESULT(PGZ_V_VW)
+! #######################################################
+!
+!!**** *GZ_V_VW - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Z
+!! direction for a variable placed at the
+!! V point and the result is placed at
+!! the VW vorticity point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the Z cartesian direction for a field PA placed at the
+! V point. The result is placed at the VW vorticity point.
+!
+!
+! dzm(PA)
+! PGZ_V_VW = ------
+! ____y
+! d*zz
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MYM : Shuman functions (mean operators)
+!! DZM : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 20/07/94
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the V point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGZ_V_VW ! result VW point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GZ_V_VW
+! ---------------------
+!
+PGZ_V_VW(:,:,:)= DZM(PA) / MYM(PDZZ)
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GZ_V_VW
diff --git a/src/mesonh/aux/gradient_w.f90 b/src/mesonh/aux/gradient_w.f90
new file mode 100644
index 000000000..1ef8f6916
--- /dev/null
+++ b/src/mesonh/aux/gradient_w.f90
@@ -0,0 +1,297 @@
+!MNH_LIC Copyright 1994-2020 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_GRADIENT_W
+! ######################
+!
+INTERFACE
+!
+!
+FUNCTION GZ_W_M(PA,PDZZ) RESULT(PGZ_W_M)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the W point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGZ_W_M ! result mass point
+!
+END FUNCTION GZ_W_M
+!
+FUNCTION GX_W_UW(PA,PDXX,PDZZ,PDZX) RESULT(PGX_W_UW)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the W point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_W_UW ! result UW point
+!
+END FUNCTION GX_W_UW
+!
+!
+FUNCTION GY_W_VW(PA,PDYY,PDZZ,PDZY) RESULT(PGY_W_VW)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the W point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_W_VW ! result VW point
+!
+END FUNCTION GY_W_VW
+!
+!
+END INTERFACE
+!
+END MODULE MODI_GRADIENT_W
+!
+!
+!
+! #######################################################
+ FUNCTION GZ_W_M(PA,PDZZ) RESULT(PGZ_W_M)
+! #######################################################
+!
+!!**** *GZ_W_M* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Z
+!! direction for a variable placed at the
+!! W point and the result is placed at
+!! the mass point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the Z cartesian direction for a field PA placed at the
+! W point. The result is placed at the mass point.
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MZF : Shuman functions (mean operators)
+!! DZF : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 19/07/94
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the W point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGZ_W_M ! result mass point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GZ_W_M
+! --------------------
+!
+PGZ_W_M(:,:,:)= DZF(PA(:,:,:))/(MZF(PDZZ(:,:,:)))
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GZ_W_M
+!
+!
+! #########################################################
+ FUNCTION GX_W_UW(PA,PDXX,PDZZ,PDZX) RESULT(PGX_W_UW)
+! #########################################################
+!
+!!**** *GX_W_UW* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian X
+!! direction for a variable placed at the
+!! V point and the result is placed at
+!! the UW vorticity point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the X cartesian direction for a field PA placed at the
+! W point. The result is placed at the UW vorticity point.
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MXM,MZM,MZF : Shuman functions (mean operators)
+!! DXM,DZM : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 20/07/94
+!! 18/10/00 (V.Masson) add LFLAT switch
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+USE MODD_CONF
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the W point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_W_UW ! result UW point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GX_W_UW
+! ---------------------
+!
+IF (.NOT. LFLAT) THEN
+ PGX_W_UW(:,:,:)= DXM(PA(:,:,:))/(MZM(PDXX(:,:,:))) &
+ -DZM(MXM(MZF(PA(:,:,:))))*PDZX(:,:,:) &
+ /( MZM(PDXX(:,:,:))*MXM(PDZZ(:,:,:)) )
+ELSE
+ PGX_W_UW(:,:,:)= DXM(PA(:,:,:))/(MZM(PDXX(:,:,:)))
+END IF
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GX_W_UW
+!
+!
+! #########################################################
+ FUNCTION GY_W_VW(PA,PDYY,PDZZ,PDZY) RESULT(PGY_W_VW)
+! #########################################################
+!
+!!**** *GY_W_VW* - Cartesian Gradient operator:
+!! computes the gradient in the cartesian Y
+!! direction for a variable placed at the
+!! W point and the result is placed at
+!! the VW vorticity point.
+!! PURPOSE
+!! -------
+! The purpose of this function is to compute the discrete gradient
+! along the Y cartesian direction for a field PA placed at the
+! W point. The result is placed at the VW vorticity point.
+!
+!!** METHOD
+!! ------
+!! The Chain rule of differencing is applied to variables expressed
+!! in the Gal-Chen & Somerville coordinates to obtain the gradient in
+!! the cartesian system
+!!
+!! EXTERNAL
+!! --------
+!! MYM,MZM,MZF : Shuman functions (mean operators)
+!! DYM,DZM : Shuman functions (finite difference operators)
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation of Meso-NH (GRAD_CAR operators)
+!! A Turbulence scheme for the Meso-NH model (Chapter 6)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart *INM and Meteo-France*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 20/07/94
+!! 18/10/00 (V.Masson) add LFLAT switch
+!-------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+!
+USE MODI_SHUMAN
+USE MODD_CONF
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at the W point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzx
+!
+REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_W_VW ! result VW point
+!
+!
+!* 0.2 declaration of local variables
+!
+! NONE
+!
+!----------------------------------------------------------------------------
+!
+!* 1. DEFINITION of GY_W_VW
+! ---------------------
+!
+IF (.NOT. LFLAT) THEN
+ PGY_W_VW(:,:,:)= DYM(PA(:,:,:))/(MZM(PDYY(:,:,:))) &
+ -DZM(MYM(MZF(PA(:,:,:))))*PDZY(:,:,:) &
+ /( MZM(PDYY(:,:,:))*MYM(PDZZ(:,:,:)) )
+ELSE
+ PGY_W_VW(:,:,:)= DYM(PA(:,:,:))/(MZM(PDYY(:,:,:)))
+END IF
+!
+!----------------------------------------------------------------------------
+!
+END FUNCTION GY_W_VW
diff --git a/src/mesonh/turb/mode_turb_ver_dyn_flux.f90 b/src/mesonh/turb/mode_turb_ver_dyn_flux.f90
new file mode 100644
index 000000000..4473628de
--- /dev/null
+++ b/src/mesonh/turb/mode_turb_ver_dyn_flux.f90
@@ -0,0 +1,778 @@
+!MNH_LIC Copyright 1994-2021 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 MODE_TURB_VER_DYN_FLUX
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE TURB_VER_DYN_FLUX(KKA,KKU,KKL, &
+ OTURB_FLX,KRR, &
+ HTURBDIM,PIMPL,PEXPL, &
+ PTSTEP, &
+ TPFILE, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ, &
+ PCOSSLOPE,PSINSLOPE, &
+ PRHODJ, &
+ PCDUEFF,PTAU11M,PTAU12M,PTAU33M, &
+ PTHLM,PRM,PSVM,PUM,PVM,PWM,PUSLOPEM,PVSLOPEM, &
+ PTKEM,PLM,MFMOIST,PWU,PWV, &
+ PRUS,PRVS,PRWS, &
+ PDP,PTP )
+! ###############################################################
+!
+!
+!!**** *TURB_VER_DYN_FLUX* -compute the source terms due to the vertical turbulent
+!! fluxes.
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute the vertical turbulent
+! fluxes of the evolutive variables and give back the source
+! terms to the main program. In the case of large horizontal meshes,
+! the divergence of these vertical turbulent fluxes represent the whole
+! effect of the turbulence but when the three-dimensionnal version of
+! the turbulence scheme is activated (CTURBDIM="3DIM"), these divergences
+! are completed in the next routine TURB_HOR.
+! An arbitrary degree of implicitness has been implemented for the
+! temporal treatment of these diffusion terms.
+! The vertical boundary conditions are as follows:
+! * at the bottom, the surface fluxes are prescribed at the same
+! as the other turbulent fluxes
+! * at the top, the turbulent fluxes are set to 0.
+! It should be noted that the condensation has been implicitely included
+! in this turbulence scheme by using conservative variables and computing
+! the subgrid variance of a statistical variable s indicating the presence
+! or not of condensation in a given mesh.
+!
+!!** METHOD
+!! ------
+!! 1D type calculations are made;
+!! The vertical turbulent fluxes are computed in an off-centered
+!! implicit scheme (a Crank-Nicholson type with coefficients different
+!! than 0.5), which allows to vary the degree of implicitness of the
+!! formulation.
+!! The different prognostic variables are treated one by one.
+!! The contributions of each turbulent fluxes are cumulated into the
+!! tendency PRvarS, and into the dynamic and thermal production of
+!! TKE if necessary.
+!!
+!! In section 2 and 3, the thermodynamical fields are considered.
+!! Only the turbulent fluxes of the conservative variables
+!! (Thetal and Rnp stored in PRx(:,:,:,1)) are computed.
+!! Note that the turbulent fluxes at the vertical
+!! boundaries are given either by the soil scheme for the surface one
+!! ( at the same instant as the others fluxes) and equal to 0 at the
+!! top of the model. The thermal production is computed by vertically
+!! averaging the turbulent flux and multiply this flux at the mass point by
+!! a function ETHETA or EMOIST, which preform the transformation from the
+!! conservative variables to the virtual potential temperature.
+!!
+!! In section 4, the variance of the statistical variable
+!! s indicating presence or not of condensation, is determined in function
+!! of the turbulent moments of the conservative variables and its
+!! squarred root is stored in PSIGS. This information will be completed in
+!! the horizontal turbulence if the turbulence dimensionality is not
+!! equal to "1DIM".
+!!
+!! In section 5, the x component of the stress tensor is computed.
+!! The surface flux <u'w'> is computed from the value of the surface
+!! fluxes computed in axes linked to the orography ( i", j" , k"):
+!! i" is parallel to the surface and in the direction of the maximum
+!! slope
+!! j" is also parallel to the surface and in the normal direction of
+!! the maximum slope
+!! k" is the normal to the surface
+!! In order to prevent numerical instability, the implicit scheme has
+!! been extended to the surface flux regarding to its dependence in
+!! function of U. The dependence in function of the other components
+!! introduced by the different rotations is only explicit.
+!! The turbulent fluxes are used to compute the dynamic production of
+!! TKE. For the last TKE level ( located at PDZZ(:,:,IKB)/2 from the
+!! ground), an harmonic extrapolation from the dynamic production at
+!! PDZZ(:,:,IKB) is used to avoid an evaluation of the gradient of U
+!! in the surface layer.
+!!
+!! In section 6, the same steps are repeated but for the y direction
+!! and in section 7, a diagnostic computation of the W variance is
+!! performed.
+!!
+!! In section 8, the turbulent fluxes for the scalar variables are
+!! computed by the same way as the conservative thermodynamical variables
+!!
+!!
+!! EXTERNAL
+!! --------
+!! GX_U_M, GY_V_M, GZ_W_M : cartesian gradient operators
+!! GX_U_UW,GY_V_VW (X,Y,Z) represent the direction of the gradient
+!! _(M,U,...)_ represent the localization of the
+!! field to be derivated
+!! _(M,UW,...) represent the localization of the
+!! field derivated
+!!
+!!
+!! MXM,MXF,MYM,MYF,MZM,MZF
+!! : Shuman functions (mean operators)
+!! DXF,DYF,DZF,DZM
+!! : Shuman functions (difference operators)
+!!
+!! SUBROUTINE TRIDIAG : to compute the split implicit evolution
+!! of a variable located at a mass point
+!!
+!! SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
+!! of a variable located at a wind point
+!!
+!! FUNCTIONs ETHETA and EMOIST :
+!! allows to compute:
+!! - the coefficients for the turbulent correlation between
+!! any variable and the virtual potential temperature, of its
+!! correlations with the conservative potential temperature and
+!! the humidity conservative variable:
+!! ------- ------- -------
+!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
+!!
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : contains physical constants
+!!
+!! XG : gravity constant
+!!
+!! Module MODD_CTURB: contains the set of constants for
+!! the turbulence scheme
+!!
+!! XCMFS,XCMFB : cts for the momentum flux
+!! XCSHF : ct for the sensible heat flux
+!! XCHF : ct for the moisture flux
+!! XCTV,XCHV : cts for the T and moisture variances
+!!
+!! Module MODD_PARAMETERS
+!!
+!! JPVEXT_TURB : number of vertical external points
+!! JPHEXT : number of horizontal external points
+!!
+!!
+!! REFERENCE
+!! ---------
+!! Book 1 of documentation (Chapter: Turbulence)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original August 19, 1994
+!! Modifications: February 14, 1995 (J.Cuxart and J.Stein)
+!! Doctorization and Optimization
+!! Modifications: March 21, 1995 (J.M. Carriere)
+!! Introduction of cloud water
+!! Modifications: June 14, 1995 (J.Cuxart and J. Stein)
+!! Phi3 and Psi3 at w-point + bug in the all
+!! or nothing condens.
+!! Modifications: Sept 15, 1995 (J.Cuxart and J. Stein)
+!! Change the DP computation at the ground
+!! Modifications: October 10, 1995 (J.Cuxart and J. Stein)
+!! Psi for scal var and LES tools
+!! Modifications: November 10, 1995 (J. Stein)
+!! change the surface relations
+!! Modifications: February 20, 1995 (J. Stein) optimization
+!! Modifications: May 21, 1996 (J. Stein)
+!! bug in the vertical flux of the V wind
+!! component for explicit computation
+!! Modifications: May 21, 1996 (N. wood)
+!! modify the computation of the vertical
+!! part or the surface tangential flux
+!! Modifications: May 21, 1996 (P. Jabouille)
+!! same modification in the Y direction
+!!
+!! Modifications: Sept 17, 1996 (J. Stein) change the moist case by using
+!! Pi instead of Piref + use Atheta and Amoist
+!!
+!! Modifications: Nov 24, 1997 (V. Masson) removes the DO loops
+!! Modifications: Mar 31, 1998 (V. Masson) splits the routine TURB_VER_DYN_FLUX
+!! Modifications: Oct 18, 2000 (J. Stein) Bug in some computations for IKB level
+!! Modifications: Oct 18, 2000 (V. Masson) LES computations + LFLAT switch
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
+!! change of YCOMMENT
+!! 2012-02 Y. Seity, add possibility to run with reversed vertical levels
+!! Modifications July 2015 (Wim de Rooy) LHARATU switch
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Q. Rodier 17/01/2019 : cleaning : remove cyclic conditions on DP and ZA
+!! JL Redelsperger 03/2021 : Add Ocean & O-A Autocoupling LES Cases
+!!--------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
+USE MODD_CONF
+USE MODD_CST
+USE MODD_CTURB
+USE MODD_DYN_n, ONLY: LOCEAN
+USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_LES
+USE MODD_NSV
+USE MODD_OCEANH
+USE MODD_PARAMETERS
+USE MODD_REF, ONLY : LCOUPLES
+USE MODD_TURB_n
+!
+!
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_GRADIENT_M
+USE MODI_SECOND_MNH
+USE MODI_SHUMAN , ONLY: MZM, MZF, MXM, MXF, MYM, MYF,&
+ & DZM, DXF, DXM, DYF, DYM
+USE MODI_TRIDIAG
+USE MODI_TRIDIAG_WIND
+USE MODI_LES_MEAN_SUBGRID
+!
+USE MODE_IO_FIELD_WRITE, only: IO_Field_write
+USE MODE_ll
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of arguments
+!
+!
+!
+INTEGER, INTENT(IN) :: KKA !near ground array index
+INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
+INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
+LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
+ ! turbulent fluxes in the syncronous FM-file
+INTEGER, INTENT(IN) :: KRR ! number of moist var.
+CHARACTER(len=4), INTENT(IN) :: HTURBDIM ! dimensionality of the
+ ! turbulence scheme
+REAL, INTENT(IN) :: PIMPL, PEXPL ! Coef. for temporal disc.
+REAL, INTENT(IN) :: PTSTEP ! Double Time Step
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSZW ! Director Cosinus of the
+ ! normal to the ground surface
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! altitude of flux points
+REAL, DIMENSION(:,:), INTENT(IN) :: PCOSSLOPE ! cosinus of the angle
+ ! between i and the slope vector
+REAL, DIMENSION(:,:), INTENT(IN) :: PSINSLOPE ! sinus of the angle
+ ! between i and the slope vector
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry density * grid volum
+REAL, DIMENSION(:,:,:), INTENT(IN) :: MFMOIST ! moist mass flux dual scheme
+
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PCDUEFF ! Cd * || u || at time t
+REAL, DIMENSION(:,:), INTENT(IN) :: PTAU11M ! <uu> in the axes linked
+ ! to the maximum slope direction and the surface normal and the binormal
+ ! at time t - dt
+REAL, DIMENSION(:,:), INTENT(IN) :: PTAU12M ! <uv> in the same axes
+REAL, DIMENSION(:,:), INTENT(IN) :: PTAU33M ! <ww> in the same axes
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUM,PVM,PWM, PTHLM
+ ! Wind at t-Delta t
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM
+REAL, DIMENSION(:,:), INTENT(IN) :: PUSLOPEM ! wind component along the
+ ! maximum slope direction
+REAL, DIMENSION(:,:), INTENT(IN) :: PVSLOPEM ! wind component along the
+ ! direction normal to the maximum slope one
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at time t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PWU ! momentum flux u'w'
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PWV ! momentum flux v'w'
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS, PRWS
+ ! cumulated sources for the prognostic variables
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDP ! Dynamic TKE production term
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTP ! Thermal TKE production term
+!
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+!
+REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2)) :: ZDIRSINZW ! sinus of the angle
+ ! between the normal and the vertical at the surface
+REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),1):: ZCOEFS ! coeff. for the
+ ! implicit scheme for the wind at the surface
+REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),SIZE(PUM,3)) :: &
+ ZA, & ! under diagonal elements of the tri-diagonal matrix involved
+ ! in the temporal implicit scheme (also used to store coefficient
+ ! J in Section 5)
+ ZRES, & ! guess of the treated variable at t+ deltat when the turbu-
+ ! lence is the only source of evolution added to the ones
+ ! considered in ZSOURCE
+ ZFLXZ, & ! vertical flux of the treated variable
+ ZSOURCE, & ! source of evolution for the treated variable
+ ZKEFF ! effectif diffusion coeff = LT * SQRT( TKE )
+INTEGER :: IRESP ! Return code of FM routines
+INTEGER :: IGRID ! C-grid indicator in LFIFM file
+INTEGER :: ILENCH ! Length of comment string in LFIFM file
+INTEGER :: IIB,IIE, & ! I index values for the Beginning and End
+ IJB,IJE, & ! mass points of the domain in the 3 direct.
+ IKB,IKE !
+INTEGER :: IKT ! array size in k direction
+INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
+INTEGER :: JSV ! scalar loop counter
+CHARACTER (LEN=100) :: YCOMMENT ! comment string in LFIFM file
+CHARACTER (LEN=16) :: YRECFM ! Name of the desired field in LFIFM file
+REAL, DIMENSION(SIZE(PDZZ,1),SIZE(PDZZ,2),1) :: ZCOEFFLXU, &
+ ZCOEFFLXV, ZUSLOPEM, ZVSLOPEM
+ ! coefficients for the surface flux
+ ! evaluation and copy of PUSLOPEM and
+ ! PVSLOPEM in local 3D arrays
+INTEGER :: IIU,IJU ! size of array in x,y,z directions
+!
+REAL :: ZTIME1, ZTIME2, ZCMFS
+TYPE(TFIELDDATA) :: TZFIELD
+!----------------------------------------------------------------------------
+!
+!* 1. PRELIMINARIES
+! -------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB_VER_DYN_FLUX',0,ZHOOK_HANDLE)
+IIU=SIZE(PUM,1)
+IIE=IIU-JPHEXT
+IIB=1+JPHEXT
+IJU=SIZE(PUM,2)
+IJE=IJU-JPHEXT
+IJB=1+JPHEXT
+IKB=KKA+JPVEXT_TURB*KKL
+IKE=KKU-JPVEXT_TURB*KKL
+IKT=SIZE(PUM,3)
+IKTB=1+JPVEXT_TURB
+IKTE=IKT-JPVEXT_TURB
+
+
+!
+ZSOURCE = 0.
+ZFLXZ = 0.
+ZCMFS = XCMFS
+IF (LHARAT)THEN
+ ZCMFS=1.
+ENDIF
+!
+ZDIRSINZW(:,:) = SQRT(1.-PDIRCOSZW(:,:)**2)
+! compute the coefficients for the uncentred gradient computation near the
+! ground
+!
+! With LHARATU length scale and TKE are at half levels so remove MZM
+!
+IF (LHARAT) THEN
+ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
+ELSE
+ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
+ENDIF
+
+!
+ZUSLOPEM(:,:,1)=PUSLOPEM(:,:)
+ZVSLOPEM(:,:,1)=PVSLOPEM(:,:)
+!
+!----------------------------------------------------------------------------
+!
+!
+!* 5. SOURCES OF U,W WIND COMPONENTS AND PARTIAL DYNAMIC PRODUCTION
+! -------------------------------------------------------------
+!
+!* 5.1 Source of U wind component
+!
+! Preparation of the arguments for TRIDIAG_WIND
+!
+ZA(:,:,:) = -PTSTEP * ZCMFS * &
+ MXM( ZKEFF ) * MXM(MZM(PRHODJ, KKA, KKU, KKL)) / &
+ MXM( PDZZ )**2
+!
+IF (CPROGRAM/='AROME ') ZA(1,:,:)=ZA(IIE,:,:)
+!
+! Compute the source of U wind component
+!
+! compute the coefficient between the vertical flux and the 2 components of the
+! wind following the slope
+ZCOEFFLXU(:,:,1) = PCDUEFF(:,:) * (PDIRCOSZW(:,:)**2 - ZDIRSINZW(:,:)**2) &
+ * PCOSSLOPE(:,:)
+ZCOEFFLXV(:,:,1) = PCDUEFF(:,:) * PDIRCOSZW(:,:) * PSINSLOPE(:,:)
+
+! prepare the implicit scheme coefficients for the surface flux
+ZCOEFS(:,:,1)= ZCOEFFLXU(:,:,1) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:) &
+ +ZCOEFFLXV(:,:,1) * PSINSLOPE(:,:)
+!
+! average this flux to be located at the U,W vorticity point
+ZCOEFS(:,:,1:1)=MXM(ZCOEFS(:,:,1:1) / PDZZ(:,:,IKB:IKB) )
+!
+! compute the explicit tangential flux at the W point
+ZSOURCE(:,:,IKB) = &
+ PTAU11M(:,:) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
+ -PTAU12M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) &
+ -PTAU33M(:,:) * PCOSSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
+!
+! add the vertical part or the surface flux at the U,W vorticity point
+
+ZSOURCE(:,:,IKB:IKB) = &
+ ( MXM( ZSOURCE(:,:,IKB:IKB) / PDZZ(:,:,IKB:IKB) ) &
+ + MXM( ZCOEFFLXU(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
+ *ZUSLOPEM(:,:,1:1) &
+ -ZCOEFFLXV(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
+ *ZVSLOPEM(:,:,1:1) ) &
+ - ZCOEFS(:,:,1:1) * PUM(:,:,IKB:IKB) * PIMPL &
+ ) * 0.5 * ( 1. + MXM(PRHODJ(:,:,KKA:KKA)) / MXM(PRHODJ(:,:,IKB:IKB)) )
+!
+ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
+ZSOURCE(:,:,IKE) = 0.
+!
+! Obtention of the splitted U at t+ deltat
+!
+CALL TRIDIAG_WIND(KKA,KKU,KKL,PUM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL, &
+ MXM(PRHODJ),ZSOURCE,ZRES)
+!
+! Compute the equivalent tendency for the U wind component
+!
+PRUS(:,:,:)=PRUS(:,:,:)+MXM(PRHODJ(:,:,:))*(ZRES(:,:,:)-PUM(:,:,:))/PTSTEP
+!
+!
+!* 5.2 Partial Dynamic Production
+!
+! vertical flux of the U wind component
+!
+ZFLXZ(:,:,:) = -ZCMFS * MXM(ZKEFF) * &
+ DZM(PIMPL*ZRES + PEXPL*PUM, KKA, KKU, KKL) / MXM(PDZZ)
+!
+! surface flux
+ZFLXZ(:,:,IKB:IKB) = MXM(PDZZ(:,:,IKB:IKB)) * &
+ ( ZSOURCE(:,:,IKB:IKB) &
+ +ZCOEFS(:,:,1:1) * ZRES(:,:,IKB:IKB) * PIMPL &
+ ) / 0.5 / ( 1. + MXM(PRHODJ(:,:,KKA:KKA)) / MXM(PRHODJ(:,:,IKB:IKB)) )
+!
+ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+
+!
+IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+ ! stores the U wind component vertical flux
+ TZFIELD%CMNHNAME = 'UW_VFLX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'UW_VFLX'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'U wind component vertical flux'
+ TZFIELD%NGRID = 4
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
+END IF
+!
+! first part of total momentum flux
+!
+PWU(:,:,:) = ZFLXZ(:,:,:)
+!
+! Contribution to the dynamic production of TKE
+! compute the dynamic production at the mass point
+!
+PDP(:,:,:) = - MZF(MXF(ZFLXZ * GZ_U_UW(PUM,PDZZ, KKA, KKU, KKL)), KKA, KKU, KKL)
+!
+! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
+PDP(:,:,IKB:IKB) = - MXF ( &
+ ZFLXZ(:,:,IKB+KKL:IKB+KKL) * (PUM(:,:,IKB+KKL:IKB+KKL)-PUM(:,:,IKB:IKB)) &
+ / MXM(PDZZ(:,:,IKB+KKL:IKB+KKL)) &
+ )
+!
+! Storage in the LES configuration
+!
+IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID(MZF(MXF(ZFLXZ), KKA, KKU, KKL), X_LES_SUBGRID_WU )
+ CALL LES_MEAN_SUBGRID(MZF(MXF(GZ_U_UW(PUM,PDZZ, KKA, KKU, KKL) &
+ & *ZFLXZ), KKA, KKU, KKL), X_LES_RES_ddxa_U_SBG_UaU )
+ CALL LES_MEAN_SUBGRID( ZCMFS * ZKEFF, X_LES_SUBGRID_Km )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+!
+!* 5.3 Source of W wind component
+!
+!
+IF(HTURBDIM=='3DIM') THEN
+ ! Compute the source for the W wind component
+ ZFLXZ(:,:,KKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation
+ ! used to compute the W source at the ground
+ !
+ IF (.NOT. LFLAT) THEN
+ PRWS(:,:,:)= PRWS &
+ -DXF( MZM(MXM(PRHODJ) /PDXX, KKA, KKU, KKL) * ZFLXZ ) &
+ +DZM(PRHODJ / MZF(PDZZ, KKA, KKU, KKL) * &
+ MXF(MZF(ZFLXZ*PDZX, KKA, KKU, KKL) / PDXX ), &
+ KKA, KKU, KKL)
+ ELSE
+ PRWS(:,:,:)= PRWS -DXF(MZM(MXM(PRHODJ) /PDXX, KKA, KKU, KKL) * ZFLXZ )
+ END IF
+ !
+ ! Complete the Dynamical production with the W wind component
+ !
+ ZA(:,:,:)=-MZF(MXF(ZFLXZ * GX_W_UW(PWM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)), KKA, KKU, KKL)
+ !
+ !
+ ! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
+ ZA(:,:,IKB:IKB) = - MXF ( &
+ ZFLXZ(:,:,IKB+KKL:IKB+KKL) * &
+ ( DXM( PWM(:,:,IKB+KKL:IKB+KKL) ) &
+ -MXM( (PWM(:,:,IKB+2*KKL:IKB+2*KKL )-PWM(:,:,IKB+KKL:IKB+KKL)) &
+ /(PDZZ(:,:,IKB+2*KKL:IKB+2*KKL)+PDZZ(:,:,IKB+KKL:IKB+KKL)) &
+ +(PWM(:,:,IKB+KKL:IKB+KKL)-PWM(:,:,IKB:IKB )) &
+ /(PDZZ(:,:,IKB+KKL:IKB+KKL)+PDZZ(:,:,IKB:IKB )) &
+ ) &
+ * PDZX(:,:,IKB+KKL:IKB+KKL) &
+ ) / (0.5*(PDXX(:,:,IKB+KKL:IKB+KKL)+PDXX(:,:,IKB:IKB))) &
+ )
+ !
+ PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
+ !
+ ! Storage in the LES configuration
+ !
+ IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID(MZF(MXF(GX_W_UW(PWM,PDXX,&
+ PDZZ,PDZX, KKA, KKU, KKL)*ZFLXZ), KKA, KKU, KKL), X_LES_RES_ddxa_W_SBG_UaW )
+ CALL LES_MEAN_SUBGRID(MXF(GX_M_U(PTHLM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)&
+ * MZF(ZFLXZ, KKA, KKU, KKL)), X_LES_RES_ddxa_Thl_SBG_UaW )
+ IF (KRR>=1) THEN
+ CALL LES_MEAN_SUBGRID(MXF(GX_U_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX, KKA, KKU, KKL)&
+ *MZF(ZFLXZ, KKA, KKU, KKL)),X_LES_RES_ddxa_Rt_SBG_UaW )
+ END IF
+ DO JSV=1,NSV
+ CALL LES_MEAN_SUBGRID( MXF(GX_U_M(PSVM(:,:,:,JSV),PDXX,PDZZ,&
+ PDZX, KKA, KKU, KKL)*MZF(ZFLXZ, KKA, KKU, KKL)),X_LES_RES_ddxa_Sv_SBG_UaW(:,:,:,JSV) )
+ END DO
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+END IF
+!
+!----------------------------------------------------------------------------
+!
+!
+!* 6. SOURCES OF V,W WIND COMPONENTS AND COMPLETE 1D DYNAMIC PRODUCTION
+! -----------------------------------------------------------------
+!
+!* 6.1 Source of V wind component
+!
+! Preparation of the arguments for TRIDIAG_WIND
+!!
+ZA(:,:,:) = - PTSTEP * ZCMFS * &
+ MYM( ZKEFF ) * MYM(MZM(PRHODJ, KKA, KKU, KKL)) / &
+ MYM( PDZZ )**2
+!
+!
+IF(CPROGRAM/='AROME ') ZA(:,1,:)=ZA(:,IJE,:)
+!
+! Compute the source of V wind component
+! compute the coefficient between the vertical flux and the 2 components of the
+! wind following the slope
+ZCOEFFLXU(:,:,1) = PCDUEFF(:,:) * (PDIRCOSZW(:,:)**2 - ZDIRSINZW(:,:)**2) &
+ * PSINSLOPE(:,:)
+ZCOEFFLXV(:,:,1) = PCDUEFF(:,:) * PDIRCOSZW(:,:) * PCOSSLOPE(:,:)
+
+! prepare the implicit scheme coefficients for the surface flux
+ZCOEFS(:,:,1)= ZCOEFFLXU(:,:,1) * PSINSLOPE(:,:) * PDIRCOSZW(:,:) &
+ +ZCOEFFLXV(:,:,1) * PCOSSLOPE(:,:)
+!
+! average this flux to be located at the V,W vorticity point
+ZCOEFS(:,:,1:1)=MYM(ZCOEFS(:,:,1:1) / PDZZ(:,:,IKB:IKB) )
+!
+! compute the explicit tangential flux at the W point
+ZSOURCE(:,:,IKB) = &
+ PTAU11M(:,:) * PSINSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
+ +PTAU12M(:,:) * PCOSSLOPE(:,:) * ZDIRSINZW(:,:) &
+ -PTAU33M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
+!
+! add the vertical part or the surface flux at the V,W vorticity point
+ZSOURCE(:,:,IKB:IKB) = &
+ ( MYM( ZSOURCE(:,:,IKB:IKB) / PDZZ(:,:,IKB:IKB) ) &
+ + MYM( ZCOEFFLXU(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
+ *ZUSLOPEM(:,:,1:1) &
+ +ZCOEFFLXV(:,:,1:1) / PDZZ(:,:,IKB:IKB) &
+ *ZVSLOPEM(:,:,1:1) ) &
+ - ZCOEFS(:,:,1:1) * PVM(:,:,IKB:IKB) * PIMPL &
+ ) * 0.5 * ( 1. + MYM(PRHODJ(:,:,KKA:KKA)) / MYM(PRHODJ(:,:,IKB:IKB)) )
+!
+ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
+ZSOURCE(:,:,IKE) = 0.
+!
+! Obtention of the splitted V at t+ deltat
+CALL TRIDIAG_WIND(KKA,KKU,KKL,PVM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL, &
+ MYM(PRHODJ),ZSOURCE,ZRES)
+!
+! Compute the equivalent tendency for the V wind component
+!
+PRVS(:,:,:)=PRVS(:,:,:)+MYM(PRHODJ(:,:,:))*(ZRES(:,:,:)-PVM(:,:,:))/PTSTEP
+!
+!
+!* 6.2 Complete 1D dynamic Production
+!
+! vertical flux of the V wind component
+!
+ZFLXZ(:,:,:) = -ZCMFS * MYM(ZKEFF) * &
+ DZM(PIMPL*ZRES + PEXPL*PVM, KKA, KKU, KKL) / MYM(PDZZ)
+!
+ZFLXZ(:,:,IKB:IKB) = MYM(PDZZ(:,:,IKB:IKB)) * &
+ ( ZSOURCE(:,:,IKB:IKB) &
+ +ZCOEFS(:,:,1:1) * ZRES(:,:,IKB:IKB) * PIMPL &
+ ) / 0.5 / ( 1. + MYM(PRHODJ(:,:,KKA:KKA)) / MYM(PRHODJ(:,:,IKB:IKB)) )
+!
+!
+ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+!
+IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+ ! stores the V wind component vertical flux
+ TZFIELD%CMNHNAME = 'VW_VFLX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'VW_VFLX'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'V wind component vertical flux'
+ TZFIELD%NGRID = 4
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
+END IF
+!
+! second part of total momentum flux
+!
+PWV(:,:,:) = ZFLXZ(:,:,:)
+!
+! Contribution to the dynamic production of TKE
+! compute the dynamic production contribution at the mass point
+!
+ZA(:,:,:) = - MZF(MYF(ZFLXZ * GZ_V_VW(PVM,PDZZ, KKA, KKU, KKL)), KKA, KKU, KKL)
+!
+! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
+ZA(:,:,IKB:IKB) = &
+ - MYF ( &
+ZFLXZ(:,:,IKB+KKL:IKB+KKL) * (PVM(:,:,IKB+KKL:IKB+KKL)-PVM(:,:,IKB:IKB)) &
+ / MYM(PDZZ(:,:,IKB+KKL:IKB+KKL)) &
+ )
+!
+PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
+!
+! Storage in the LES configuration
+!
+IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID(MZF(MYF(ZFLXZ), KKA, KKU, KKL), X_LES_SUBGRID_WV )
+ CALL LES_MEAN_SUBGRID(MZF(MYF(GZ_V_VW(PVM,PDZZ, KKA, KKU, KKL)*&
+ & ZFLXZ), KKA, KKU, KKL), X_LES_RES_ddxa_V_SBG_UaV )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+!
+!
+!* 6.3 Source of W wind component
+!
+IF(HTURBDIM=='3DIM') THEN
+ ! Compute the source for the W wind component
+ ZFLXZ(:,:,KKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation
+ !
+ IF (.NOT. L2D) THEN
+ IF (.NOT. LFLAT) THEN
+ PRWS(:,:,:)= PRWS(:,:,:) &
+ -DYF( MZM(MYM(PRHODJ) /PDYY, KKA, KKU, KKL) * ZFLXZ ) &
+ +DZM(PRHODJ / MZF(PDZZ, KKA, KKU, KKL) * &
+ MYF(MZF(ZFLXZ*PDZY, KKA, KKU, KKL) / PDYY ), &
+ KKA, KKU, KKL)
+ ELSE
+ PRWS(:,:,:)= PRWS(:,:,:) -DYF(MZM(MYM(PRHODJ) /PDYY, KKA, KKU, KKL) * ZFLXZ )
+ END IF
+ END IF
+ !
+ ! Complete the Dynamical production with the W wind component
+ IF (.NOT. L2D) THEN
+ ZA(:,:,:) = - MZF(MYF(ZFLXZ * GY_W_VW(PWM,PDYY,PDZZ,PDZY, KKA, KKU, KKL)), KKA, KKU, KKL)
+ !
+ ! evaluate the dynamic production at w(IKB+KKL) in PDP(IKB)
+ ZA(:,:,IKB:IKB) = - MYF ( &
+ ZFLXZ(:,:,IKB+KKL:IKB+KKL) * &
+ ( DYM( PWM(:,:,IKB+KKL:IKB+KKL) ) &
+ -MYM( (PWM(:,:,IKB+2*KKL:IKB+2*KKL)-PWM(:,:,IKB+KKL:IKB+KKL)) &
+ /(PDZZ(:,:,IKB+2*KKL:IKB+2*KKL)+PDZZ(:,:,IKB+KKL:IKB+KKL)) &
+ +(PWM(:,:,IKB+KKL:IKB+KKL)-PWM(:,:,IKB:IKB )) &
+ /(PDZZ(:,:,IKB+KKL:IKB+KKL)+PDZZ(:,:,IKB:IKB )) &
+ ) &
+ * PDZY(:,:,IKB+KKL:IKB+KKL) &
+ ) / (0.5*(PDYY(:,:,IKB+KKL:IKB+KKL)+PDYY(:,:,IKB:IKB))) &
+ )
+ !
+ PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
+ !
+ END IF
+ !
+ ! Storage in the LES configuration
+ !
+ IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID(MZF(MYF(GY_W_VW(PWM,PDYY,&
+ &PDZZ,PDZY, KKA, KKU, KKL)*ZFLXZ), KKA, KKU, KKL), &
+ &X_LES_RES_ddxa_W_SBG_UaW , .TRUE. )
+ CALL LES_MEAN_SUBGRID(MYF(GY_M_V(PTHLM,PDYY,PDZZ,PDZY, KKA, KKU, KKL)*&
+ &MZF(ZFLXZ, KKA, KKU, KKL)), &
+ &X_LES_RES_ddxa_Thl_SBG_UaW , .TRUE. )
+ IF (KRR>=1) THEN
+ CALL LES_MEAN_SUBGRID(MYF(GY_V_M(PRM(:,:,:,1),PDYY,PDZZ,&
+ &PDZY, KKA, KKU, KKL)*MZF(ZFLXZ, KKA, KKU, KKL)),&
+ &X_LES_RES_ddxa_Rt_SBG_UaW , .TRUE. )
+ END IF
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+ !
+END IF
+!
+! complete the dynamic production at the marginal points
+IF (CPROGRAM/='AROME ') THEN
+ PDP(:,:,KKA)= -999.
+ PDP(:,:,KKU)= -999.
+ PDP(:,1,:)= PDP(:,IJE,:)
+ PDP(:,IJE+1,:)= PDP(:,IJB,:)
+ PDP(1,:,:)= PDP(IIE,:,:)
+ PDP(IIE+1,:,:)= PDP(IIB,:,:)
+END IF
+!
+!----------------------------------------------------------------------------
+!
+!* 7. DIAGNOSTIC COMPUTATION OF THE 1D <W W> VARIANCE
+! -----------------------------------------------
+!
+IF ( OTURB_FLX .AND. TPFILE%LOPENED .AND. HTURBDIM == '1DIM') THEN
+ ZFLXZ(:,:,:)= (2./3.) * PTKEM(:,:,:) &
+ -ZCMFS*PLM(:,:,:)*SQRT(PTKEM(:,:,:))*GZ_W_M(PWM,PDZZ, KKA, KKU, KKL)
+ ! to be tested &
+ ! +XCMFB*(4./3.)*PLM(:,:,:)/SQRT(PTKEM(:,:,:))*PTP(:,:,:)
+ ! stores the W variance
+ TZFIELD%CMNHNAME = 'W_VVAR'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'W_VVAR'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_W_VVAR'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
+END IF
+!
+!----------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('TURB_VER_DYN_FLUX',1,ZHOOK_HANDLE)
+END SUBROUTINE TURB_VER_DYN_FLUX
+END MODULE MODE_TURB_VER_DYN_FLUX
diff --git a/src/mesonh/turb/mode_turb_ver_sv_corr.f90 b/src/mesonh/turb/mode_turb_ver_sv_corr.f90
new file mode 100644
index 000000000..f140c0792
--- /dev/null
+++ b/src/mesonh/turb/mode_turb_ver_sv_corr.f90
@@ -0,0 +1,178 @@
+!MNH_LIC Copyright 2002-2020 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 MODE_TURB_VER_SV_CORR
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE TURB_VER_SV_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI, &
+ PDZZ, &
+ PTHLM,PRM,PTHVREF, &
+ PLOCPEXNM,PATHETA,PAMOIST,PSRCM,PPHI3,PPSI3, &
+ PWM,PSVM, &
+ PTKEM,PLM,PLEPS,PPSI_SV )
+! ###############################################################
+!
+!
+!!**** *TURB_VER_SV_CORR* -compute the subgrid Sv2 and SvThv terms
+!!
+!! PURPOSE
+!! -------
+!!
+!!
+!! EXTERNAL
+!! --------
+!!
+!! FUNCTIONs ETHETA and EMOIST :
+!! allows to compute:
+!! - the coefficients for the turbulent correlation between
+!! any variable and the virtual potential temperature, of its
+!! correlations with the conservative potential temperature and
+!! the humidity conservative variable:
+!! ------- ------- -------
+!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
+!!
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original October 29, 2002
+!! JP Pinty Feb 20, 2003 Add PFRAC_ICE
+!!--------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
+USE MODD_CST
+USE MODD_CTURB
+USE MODD_PARAMETERS
+USE MODD_LES
+USE MODD_CONF
+USE MODD_NSV, ONLY : NSV,NSV_LGBEG,NSV_LGEND
+USE MODD_BLOWSNOW
+!
+!
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_GRADIENT_M
+USE MODI_SHUMAN , ONLY : MZF
+USE MODI_EMOIST
+USE MODI_ETHETA
+USE MODI_LES_MEAN_SUBGRID
+!
+USE MODI_SECOND_MNH
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of arguments
+!
+!
+!
+INTEGER, INTENT(IN) :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN) :: KKL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
+INTEGER, INTENT(IN) :: KRR ! number of moist var.
+INTEGER, INTENT(IN) :: KRRL ! number of liquid var.
+INTEGER, INTENT(IN) :: KRRI ! number of ice var.
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ
+ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM ! potential temperature at t-Delta t
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! Mixing ratios at t-Delta t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! reference Thv
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLOCPEXNM ! Lv(T)/Cp/Exnref at time t-1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PATHETA ! coefficients between
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PAMOIST ! s and Thetal and Rnp
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRCM ! normalized
+ ! 2nd-order flux s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPHI3 ! Inv.Turb.Sch.for temperature
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPSI3 ! Inv.Turb.Sch.for humidity
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM ! w at time t
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! scalar var. at t-Delta t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at time t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLEPS ! dissipative length
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PPSI_SV ! Inv.Turb.Sch.for scalars
+ ! cumulated sources for the prognostic variables
+!
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+!
+REAL, DIMENSION(SIZE(PSVM,1),SIZE(PSVM,2),SIZE(PSVM,3)) :: &
+ ZA, ZFLXZ
+!
+REAL :: ZCSV !constant for the scalar flux
+!
+INTEGER :: JSV ! loop counters
+!
+REAL :: ZTIME1, ZTIME2
+!
+REAL :: ZCSVD = 1.2 ! constant for scalar variance dissipation
+REAL :: ZCTSVD = 2.4 ! constant for temperature - scalar covariance dissipation
+REAL :: ZCQSVD = 2.4 ! constant for humidity - scalar covariance dissipation
+!----------------------------------------------------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB_VER_SV_CORR',0,ZHOOK_HANDLE)
+CALL SECOND_MNH(ZTIME1)
+!
+DO JSV=1,NSV
+ !
+ IF (LNOMIXLG .AND. JSV >= NSV_LGBEG .AND. JSV<= NSV_LGEND) CYCLE
+ !
+ ! variance Sv2
+ !
+ IF (LLES_CALL) THEN
+ ! approximation: diagnosed explicitely (without implicit term)
+ ZFLXZ(:,:,:) = PPSI_SV(:,:,:,JSV)*GZ_M_W(PSVM(:,:,:,JSV),PDZZ, KKA, KKU, KKL)**2
+ ZFLXZ(:,:,:) = XCHF / ZCSVD * PLM * PLEPS * MZF(ZFLXZ(:,:,:), KKA, KKU, KKL)
+ CALL LES_MEAN_SUBGRID(-2.*ZCSVD*SQRT(PTKEM)*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Sv2(:,:,:,JSV) )
+ CALL LES_MEAN_SUBGRID(MZF(PWM, KKA, KKU, KKL)*ZFLXZ, X_LES_RES_W_SBG_Sv2(:,:,:,JSV) )
+ END IF
+ !
+ ! covariance ThvSv
+ !
+ IF (LLES_CALL) THEN
+ ! approximation: diagnosed explicitely (without implicit term)
+ ZA(:,:,:) = ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM)
+ ZFLXZ(:,:,:)= ( XCSHF * PPHI3 + XCHF * PPSI_SV(:,:,:,JSV) ) &
+ * GZ_M_W(PTHLM,PDZZ, KKA, KKU, KKL) &
+ * GZ_M_W(PSVM(:,:,:,JSV),PDZZ, KKA, KKU, KKL)
+ ZFLXZ(:,:,:)= PLM * PLEPS / (2.*ZCTSVD) * MZF(ZFLXZ, KKA, KKU, KKL)
+ CALL LES_MEAN_SUBGRID( ZA*ZFLXZ, X_LES_SUBGRID_SvThv(:,:,:,JSV) )
+ CALL LES_MEAN_SUBGRID( -XG/PTHVREF/3.*ZA*ZFLXZ, X_LES_SUBGRID_SvPz(:,:,:,JSV), .TRUE.)
+ !
+ IF (KRR>=1) THEN
+ ZA(:,:,:) = EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM)
+ ZFLXZ(:,:,:)= ( XCHF * PPSI3 + XCHF * PPSI_SV(:,:,:,JSV) ) &
+ * GZ_M_W(PRM(:,:,:,1),PDZZ, KKA, KKU, KKL) &
+ * GZ_M_W(PSVM(:,:,:,JSV),PDZZ, KKA, KKU, KKL)
+ ZFLXZ(:,:,:)= PLM * PLEPS / (2.*ZCQSVD) * MZF(ZFLXZ, KKA, KKU, KKL)
+ CALL LES_MEAN_SUBGRID( ZA*ZFLXZ, X_LES_SUBGRID_SvThv(:,:,:,JSV) , .TRUE.)
+ CALL LES_MEAN_SUBGRID( -XG/PTHVREF/3.*ZA*ZFLXZ, X_LES_SUBGRID_SvPz(:,:,:,JSV), .TRUE.)
+ END IF
+ END IF
+ !
+END DO ! end of scalar loop
+!
+CALL SECOND_MNH(ZTIME2)
+XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+!----------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('TURB_VER_SV_CORR',1,ZHOOK_HANDLE)
+END SUBROUTINE TURB_VER_SV_CORR
+END MODULE MODE_TURB_VER_SV_CORR
diff --git a/src/mesonh/turb/mode_turb_ver_sv_flux.f90 b/src/mesonh/turb/mode_turb_ver_sv_flux.f90
new file mode 100644
index 000000000..9ef220cf8
--- /dev/null
+++ b/src/mesonh/turb/mode_turb_ver_sv_flux.f90
@@ -0,0 +1,446 @@
+!MNH_LIC Copyright 1994-2021 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 MODE_TURB_VER_SV_FLUX
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE TURB_VER_SV_FLUX(KKA,KKU,KKL, &
+ OTURB_FLX,HTURBDIM, &
+ PIMPL,PEXPL, &
+ PTSTEP, &
+ TPFILE, &
+ PDZZ,PDIRCOSZW, &
+ PRHODJ,PWM, &
+ PSFSVM,PSFSVP, &
+ PSVM, &
+ PTKEM,PLM,MFMOIST,PPSI_SV, &
+ PRSVS,PWSV )
+!
+!
+!
+!!**** *TURB_VER_SV_FLUX* -compute the source terms due to the vertical turbulent
+!! fluxes.
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute the vertical turbulent
+! fluxes of the evolutive variables and give back the source
+! terms to the main program. In the case of large horizontal meshes,
+! the divergence of these vertical turbulent fluxes represent the whole
+! effect of the turbulence but when the three-dimensionnal version of
+! the turbulence scheme is activated (CTURBDIM="3DIM"), these divergences
+! are completed in the next routine TURB_HOR.
+! An arbitrary degree of implicitness has been implemented for the
+! temporal treatment of these diffusion terms.
+! The vertical boundary conditions are as follows:
+! * at the bottom, the surface fluxes are prescribed at the same
+! as the other turbulent fluxes
+! * at the top, the turbulent fluxes are set to 0.
+! It should be noted that the condensation has been implicitely included
+! in this turbulence scheme by using conservative variables and computing
+! the subgrid variance of a statistical variable s indicating the presence
+! or not of condensation in a given mesh.
+!
+!!** METHOD
+!! ------
+!! 1D type calculations are made;
+!! The vertical turbulent fluxes are computed in an off-centered
+!! implicit scheme (a Crank-Nicholson type with coefficients different
+!! than 0.5), which allows to vary the degree of implicitness of the
+!! formulation.
+!! The different prognostic variables are treated one by one.
+!! The contributions of each turbulent fluxes are cumulated into the
+!! tendency PRvarS, and into the dynamic and thermal production of
+!! TKE if necessary.
+!!
+!! In section 2 and 3, the thermodynamical fields are considered.
+!! Only the turbulent fluxes of the conservative variables
+!! (Thetal and Rnp stored in PRx(:,:,:,1)) are computed.
+!! Note that the turbulent fluxes at the vertical
+!! boundaries are given either by the soil scheme for the surface one
+!! ( at the same instant as the others fluxes) and equal to 0 at the
+!! top of the model. The thermal production is computed by vertically
+!! averaging the turbulent flux and multiply this flux at the mass point by
+!! a function ETHETA or EMOIST, which preform the transformation from the
+!! conservative variables to the virtual potential temperature.
+!!
+!! In section 4, the variance of the statistical variable
+!! s indicating presence or not of condensation, is determined in function
+!! of the turbulent moments of the conservative variables and its
+!! squarred root is stored in PSIGS. This information will be completed in
+!! the horizontal turbulence if the turbulence dimensionality is not
+!! equal to "1DIM".
+!!
+!! In section 5, the x component of the stress tensor is computed.
+!! The surface flux <u'w'> is computed from the value of the surface
+!! fluxes computed in axes linked to the orography ( i", j" , k"):
+!! i" is parallel to the surface and in the direction of the maximum
+!! slope
+!! j" is also parallel to the surface and in the normal direction of
+!! the maximum slope
+!! k" is the normal to the surface
+!! In order to prevent numerical instability, the implicit scheme has
+!! been extended to the surface flux regarding to its dependence in
+!! function of U. The dependence in function of the other components
+!! introduced by the different rotations is only explicit.
+!! The turbulent fluxes are used to compute the dynamic production of
+!! TKE. For the last TKE level ( located at PDZZ(:,:,IKB)/2 from the
+!! ground), an harmonic extrapolation from the dynamic production at
+!! PDZZ(:,:,IKB) is used to avoid an evaluation of the gradient of U
+!! in the surface layer.
+!!
+!! In section 6, the same steps are repeated but for the y direction
+!! and in section 7, a diagnostic computation of the W variance is
+!! performed.
+!!
+!! In section 8, the turbulent fluxes for the scalar variables are
+!! computed by the same way as the conservative thermodynamical variables
+!!
+!!
+!! EXTERNAL
+!! --------
+!! GX_U_M, GY_V_M, GZ_W_M : cartesian gradient operators
+!! GX_U_UW,GY_V_VW (X,Y,Z) represent the direction of the gradient
+!! _(M,U,...)_ represent the localization of the
+!! field to be derivated
+!! _(M,UW,...) represent the localization of the
+!! field derivated
+!!
+!!
+!! MXM,MXF,MYM,MYF,MZM,MZF
+!! : Shuman functions (mean operators)
+!! DXF,DYF,DZF,DZM
+!! : Shuman functions (difference operators)
+!!
+!! SUBROUTINE TRIDIAG : to compute the split implicit evolution
+!! of a variable located at a mass point
+!!
+!! SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
+!! of a variable located at a wind point
+!!
+!! FUNCTIONs ETHETA and EMOIST :
+!! allows to compute:
+!! - the coefficients for the turbulent correlation between
+!! any variable and the virtual potential temperature, of its
+!! correlations with the conservative potential temperature and
+!! the humidity conservative variable:
+!! ------- ------- -------
+!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
+!!
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : contains physical constants
+!!
+!! XG : gravity constant
+!!
+!! Module MODD_CTURB: contains the set of constants for
+!! the turbulence scheme
+!!
+!! XCMFS,XCMFB : cts for the momentum flux
+!! XCSHF : ct for the sensible heat flux
+!! XCHF : ct for the moisture flux
+!! XCTV,XCHV : cts for the T and moisture variances
+!!
+!! Module MODD_PARAMETERS
+!!
+!! JPVEXT_TURB : number of vertical external points
+!! JPHEXT : number of horizontal external points
+!!
+!!
+!! REFERENCE
+!! ---------
+!! Book 1 of documentation (Chapter: Turbulence)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original August 19, 1994
+!! Modifications: February 14, 1995 (J.Cuxart and J.Stein)
+!! Doctorization and Optimization
+!! Modifications: March 21, 1995 (J.M. Carriere)
+!! Introduction of cloud water
+!! Modifications: June 14, 1995 (J.Cuxart and J. Stein)
+!! Phi3 and Psi3 at w-point + bug in the all
+!! or nothing condens.
+!! Modifications: Sept 15, 1995 (J.Cuxart and J. Stein)
+!! Change the DP computation at the ground
+!! Modifications: October 10, 1995 (J.Cuxart and J. Stein)
+!! Psi for scal var and LES tools
+!! Modifications: November 10, 1995 (J. Stein)
+!! change the surface relations
+!! Modifications: February 20, 1995 (J. Stein) optimization
+!! Modifications: May 21, 1996 (J. Stein)
+!! bug in the vertical flux of the V wind
+!! component for explicit computation
+!! Modifications: May 21, 1996 (N. wood)
+!! modify the computation of the vertical
+!! part or the surface tangential flux
+!! Modifications: May 21, 1996 (P. Jabouille)
+!! same modification in the Y direction
+!!
+!! Modifications: Sept 17, 1996 (J. Stein) change the moist case by using
+!! Pi instead of Piref + use Atheta and Amoist
+!!
+!! Modifications: Nov 24, 1997 (V. Masson) removes the DO loops
+!! Modifications: Mar 31, 1998 (V. Masson) splits the routine TURB_VER_SV_FLUX
+!! Modifications: Dec 01, 2000 (V. Masson) conservation of scalar emission
+!! from surface in 1DIM case
+!! when slopes are present
+!! Jun 20, 2001 (J Stein) case of lagragian variables
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
+!! change of YCOMMENT
+!! Feb 2012(Y. Seity) add possibility to run with reversed
+!! vertical levels
+!! Modifications: July 2015 (Wim de Rooy) LHARAT switch
+!! Feb 2017(M. Leriche) add initialisation of ZSOURCE
+!! to avoid unknwon values outside physical domain
+!! and avoid negative values in sv tendencies
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!!--------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
+USE MODD_CST
+USE MODD_CTURB
+USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_PARAMETERS
+USE MODD_LES
+USE MODD_CONF
+USE MODD_NSV, ONLY: XSVMIN, NSV_LGBEG, NSV_LGEND
+USE MODD_BLOWSNOW
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
+!
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_GRADIENT_M
+USE MODI_SHUMAN , ONLY : DZM, MZM, MZF
+USE MODI_TRIDIAG
+USE MODI_TRIDIAG_WIND
+USE MODI_EMOIST
+USE MODI_ETHETA
+USE MODI_LES_MEAN_SUBGRID
+!
+USE MODI_SECOND_MNH
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of arguments
+!
+!
+INTEGER, INTENT(IN) :: KKA !near ground array index
+INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
+INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
+LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
+ ! turbulent fluxes in the syncronous FM-file
+CHARACTER(len=4), INTENT(IN) :: HTURBDIM ! dimensionality of the
+ ! turbulence scheme
+REAL, INTENT(IN) :: PIMPL, PEXPL ! Coef. for temporal disc.
+REAL, INTENT(IN) :: PTSTEP ! Double Time Step
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ
+ ! Metric coefficients
+REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSZW ! Director Cosinus of the
+ ! normal to the ground surface
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry density * grid volum
+REAL, DIMENSION(:,:,:), INTENT(IN) :: MFMOIST ! moist mf dual scheme
+
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSFSVM ! t - deltat
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSFSVP ! t + deltat
+!
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! scalar var. at t-Delta t
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM ! vertical wind
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at time t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PPSI_SV ! Inv.Turb.Sch.for scalars
+!
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRSVS
+ ! cumulated sources for the prognostic variables
+REAL, DIMENSION(:,:,:,:), INTENT(OUT) :: PWSV ! scalar flux
+!
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+!
+REAL, DIMENSION(SIZE(PSVM,1),SIZE(PSVM,2),SIZE(PSVM,3)) :: &
+ ZA, & ! under diagonal elements of the tri-diagonal matrix involved
+ ! in the temporal implicit scheme (also used to store coefficient
+ ! J in Section 5)
+ ZRES, & ! guess of the treated variable at t+ deltat when the turbu-
+ ! lence is the only source of evolution added to the ones
+ ! considered in ZSOURCE
+ ZFLXZ, & ! vertical flux of the treated variable
+ ZSOURCE, & ! source of evolution for the treated variable
+ ZKEFF ! effectif diffusion coeff = LT * SQRT( TKE )
+INTEGER :: IRESP ! Return code of FM routines
+INTEGER :: IGRID ! C-grid indicator in LFIFM file
+INTEGER :: ILENCH ! Length of comment string in LFIFM file
+INTEGER :: IKB,IKE ! I index values for the Beginning and End
+ ! mass points of the domain in the 3 direct.
+INTEGER :: IKT ! array size in k direction
+INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
+INTEGER :: JSV ! loop counters
+INTEGER :: JK ! loop
+INTEGER :: ISV ! number of scalar var.
+CHARACTER (LEN=100) :: YCOMMENT ! comment string in LFIFM file
+CHARACTER (LEN=16) :: YRECFM ! Name of the desired field in LFIFM file
+!
+REAL :: ZTIME1, ZTIME2
+
+REAL :: ZCSVP = 4.0 ! constant for scalar flux presso-correlation (RS81)
+REAL :: ZCSV !constant for the scalar flux
+!
+TYPE(TFIELDDATA) :: TZFIELD
+!----------------------------------------------------------------------------
+!
+!* 1. PRELIMINARIES
+! -------------
+!
+
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB_VER_SV_FLUX',0,ZHOOK_HANDLE)
+IKB=KKA+JPVEXT_TURB*KKL
+IKE=KKU-JPVEXT_TURB*KKL
+IKT=SIZE(PSVM,3)
+IKTE =IKT-JPVEXT_TURB
+IKTB =1+JPVEXT_TURB
+!
+ISV=SIZE(PSVM,4)
+!
+IF (LHARAT) THEN
+ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50.*MFMOIST(:,:,:)
+ELSE
+ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
+ENDIF
+
+!
+!----------------------------------------------------------------------------
+!
+!* 8. SOURCES OF PASSIVE SCALAR VARIABLES
+! -----------------------------------
+!
+DO JSV=1,ISV
+!
+ IF (LNOMIXLG .AND. JSV >= NSV_LGBEG .AND. JSV<= NSV_LGEND) CYCLE
+!
+! Preparation of the arguments for TRIDIAG
+IF (LHARAT) THEN
+ ZA(:,:,:) = -PTSTEP* &
+ ZKEFF * MZM(PRHODJ, KKA, KKU, KKL) / &
+ PDZZ**2
+ELSE
+ ZA(:,:,:) = -PTSTEP*XCHF*PPSI_SV(:,:,:,JSV) * &
+ ZKEFF * MZM(PRHODJ, KKA, KKU, KKL) / &
+ PDZZ**2
+ENDIF
+!
+! Compute the sources for the JSVth scalar variable
+
+!* in 3DIM case, a part of the flux goes vertically, and another goes horizontally
+! (in presence of slopes)
+!* in 1DIM case, the part of energy released in horizontal flux
+! is taken into account in the vertical part
+ IF (HTURBDIM=='3DIM') THEN
+ ZSOURCE(:,:,IKB) = (PIMPL*PSFSVP(:,:,JSV) + PEXPL*PSFSVM(:,:,JSV)) / &
+ PDZZ(:,:,IKB) * PDIRCOSZW(:,:) &
+ * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
+ ELSE
+
+ ZSOURCE(:,:,IKB) = (PIMPL*PSFSVP(:,:,JSV) + PEXPL*PSFSVM(:,:,JSV)) / &
+ PDZZ(:,:,IKB) / PDIRCOSZW(:,:) &
+ * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
+ END IF
+ ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
+ ZSOURCE(:,:,IKE) = 0.
+!
+! Obtention of the split JSV scalar variable at t+ deltat
+ CALL TRIDIAG(KKA,KKU,KKL,PSVM(:,:,:,JSV),ZA,PTSTEP,PEXPL,PIMPL,PRHODJ,ZSOURCE,ZRES)
+!
+! Compute the equivalent tendency for the JSV scalar variable
+ PRSVS(:,:,:,JSV)= PRSVS(:,:,:,JSV)+ &
+ PRHODJ(:,:,:)*(ZRES(:,:,:)-PSVM(:,:,:,JSV))/PTSTEP
+!
+ IF ( (OTURB_FLX .AND. TPFILE%LOPENED) .OR. LLES_CALL ) THEN
+ ! Diagnostic of the cartesian vertical flux
+ !
+ ZFLXZ(:,:,:) = -XCHF * PPSI_SV(:,:,:,JSV) * MZM(PLM*SQRT(PTKEM), KKA, KKU, KKL) / PDZZ * &
+ DZM(PIMPL*ZRES(:,:,:) + PEXPL*PSVM(:,:,:,JSV), KKA, KKU, KKL)
+ ! surface flux
+ !* in 3DIM case, a part of the flux goes vertically, and another goes horizontally
+ ! (in presence of slopes)
+ !* in 1DIM case, the part of energy released in horizontal flux
+ ! is taken into account in the vertical part
+ IF (HTURBDIM=='3DIM') THEN
+ ZFLXZ(:,:,IKB) = (PIMPL*PSFSVP(:,:,JSV) + PEXPL*PSFSVM(:,:,JSV)) &
+ * PDIRCOSZW(:,:)
+ ELSE
+ ZFLXZ(:,:,IKB) = (PIMPL*PSFSVP(:,:,JSV) + PEXPL*PSFSVM(:,:,JSV)) &
+ / PDIRCOSZW(:,:)
+ END IF
+ ! extrapolates the flux under the ground so that the vertical average with
+ ! the IKB flux gives the ground value
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ DO JK=IKTB+1,IKTE-1
+ PWSV(:,:,JK,JSV)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+KKL))
+ END DO
+ PWSV(:,:,IKB,JSV)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+KKL))
+ PWSV(:,:,IKE,JSV)=PWSV(:,:,IKE-KKL,JSV)
+ END IF
+ !
+ IF (OTURB_FLX .AND. TPFILE%LOPENED) THEN
+ ! stores the JSVth vertical flux
+ WRITE(TZFIELD%CMNHNAME,'("WSV_FLX_",I3.3)') JSV
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ !PW: TODO: use the correct units of the JSV variable (and multiply it by m s-1)
+ TZFIELD%CUNITS = 'SVUNIT m s-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_'//TRIM(TZFIELD%CMNHNAME)
+ TZFIELD%NGRID = 4
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ !
+ CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
+ END IF
+ !
+ ! Storage in the LES configuration
+ !
+ IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WSv(:,:,:,JSV) )
+ CALL LES_MEAN_SUBGRID(GZ_W_M(PWM,PDZZ, KKA, KKU, KKL)*MZF(ZFLXZ, KKA, KKU, KKL), &
+ X_LES_RES_ddxa_W_SBG_UaSv(:,:,:,JSV) )
+ CALL LES_MEAN_SUBGRID(MZF(GZ_M_W(PSVM(:,:,:,JSV),PDZZ, KKA, KKU, KKL)*ZFLXZ, KKA, KKU, KKL), &
+ X_LES_RES_ddxa_Sv_SBG_UaSv(:,:,:,JSV) )
+ CALL LES_MEAN_SUBGRID(-ZCSVP*SQRT(PTKEM)/PLM*MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_SvPz(:,:,:,JSV) )
+ CALL LES_MEAN_SUBGRID(MZF(PWM*ZFLXZ, KKA, KKU, KKL), X_LES_RES_W_SBG_WSv(:,:,:,JSV) )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+ !
+END DO ! end of scalar loop
+!
+!----------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('TURB_VER_SV_FLUX',1,ZHOOK_HANDLE)
+END SUBROUTINE TURB_VER_SV_FLUX
+END MODULE MODE_TURB_VER_SV_FLUX
diff --git a/src/mesonh/turb/mode_turb_ver_thermo_corr.f90 b/src/mesonh/turb/mode_turb_ver_thermo_corr.f90
new file mode 100644
index 000000000..3d420f393
--- /dev/null
+++ b/src/mesonh/turb/mode_turb_ver_thermo_corr.f90
@@ -0,0 +1,866 @@
+!MNH_LIC Copyright 1994-2021 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 MODE_TURB_VER_THERMO_CORR
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI, &
+ OTURB_FLX,HTURBDIM,HTOM, &
+ PIMPL,PEXPL, &
+ TPFILE, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW, &
+ PRHODJ,PTHVREF, &
+ PSFTHM,PSFRM,PSFTHP,PSFRP, &
+ PWM,PTHLM,PRM,PSVM, &
+ PTKEM,PLM,PLEPS, &
+ PLOCPEXNM,PATHETA,PAMOIST,PSRCM, &
+ PBETA, PSQRT_TKE, PDTH_DZ, PDR_DZ, PRED2TH3, &
+ PRED2R3, PRED2THR3, PBLL_O_E, PETHETA, &
+ PEMOIST, PREDTH1, PREDR1, PPHI3, PPSI3, PD, &
+ PFWTH,PFWR,PFTH2,PFR2,PFTHR, &
+ PTHLP,PRP,MFMOIST,PSIGS )
+! ###############################################################
+!
+!
+!!**** *TURB_VER_THERMO_FLUX* -compute the source terms due to the vertical turbulent
+!! fluxes.
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute the vertical turbulent
+! fluxes of the evolutive variables and give back the source
+! terms to the main program. In the case of large horizontal meshes,
+! the divergence of these vertical turbulent fluxes represent the whole
+! effect of the turbulence but when the three-dimensionnal version of
+! the turbulence scheme is activated (CTURBDIM="3DIM"), these divergences
+! are completed in the next routine TURB_HOR.
+! An arbitrary degree of implicitness has been implemented for the
+! temporal treatment of these diffusion terms.
+! The vertical boundary conditions are as follows:
+! * at the bottom, the surface fluxes are prescribed at the same
+! as the other turbulent fluxes
+! * at the top, the turbulent fluxes are set to 0.
+! It should be noted that the condensation has been implicitely included
+! in this turbulence scheme by using conservative variables and computing
+! the subgrid variance of a statistical variable s indicating the presence
+! or not of condensation in a given mesh.
+!
+!!** METHOD
+!! ------
+!! 1D type calculations are made;
+!! The vertical turbulent fluxes are computed in an off-centered
+!! implicit scheme (a Crank-Nicholson type with coefficients different
+!! than 0.5), which allows to vary the degree of implicitness of the
+!! formulation.
+!! The different prognostic variables are treated one by one.
+!! The contributions of each turbulent fluxes are cumulated into the
+!! tendency PRvarS, and into the dynamic and thermal production of
+!! TKE if necessary.
+!!
+!! In section 2 and 3, the thermodynamical fields are considered.
+!! Only the turbulent fluxes of the conservative variables
+!! (Thetal and Rnp stored in PRx(:,:,:,1)) are computed.
+!! Note that the turbulent fluxes at the vertical
+!! boundaries are given either by the soil scheme for the surface one
+!! ( at the same instant as the others fluxes) and equal to 0 at the
+!! top of the model. The thermal production is computed by vertically
+!! averaging the turbulent flux and multiply this flux at the mass point by
+!! a function ETHETA or EMOIST, which preform the transformation from the
+!! conservative variables to the virtual potential temperature.
+!!
+!! In section 4, the variance of the statistical variable
+!! s indicating presence or not of condensation, is determined in function
+!! of the turbulent moments of the conservative variables and its
+!! squarred root is stored in PSIGS. This information will be completed in
+!! the horizontal turbulence if the turbulence dimensionality is not
+!! equal to "1DIM".
+!!
+!! In section 5, the x component of the stress tensor is computed.
+!! The surface flux <u'w'> is computed from the value of the surface
+!! fluxes computed in axes linked to the orography ( i", j" , k"):
+!! i" is parallel to the surface and in the direction of the maximum
+!! slope
+!! j" is also parallel to the surface and in the normal direction of
+!! the maximum slope
+!! k" is the normal to the surface
+!! In order to prevent numerical instability, the implicit scheme has
+!! been extended to the surface flux regarding to its dependence in
+!! function of U. The dependence in function of the other components
+!! introduced by the different rotations is only explicit.
+!! The turbulent fluxes are used to compute the dynamic production of
+!! TKE. For the last TKE level ( located at PDZZ(:,:,IKB)/2 from the
+!! ground), an harmonic extrapolation from the dynamic production at
+!! PDZZ(:,:,IKB) is used to avoid an evaluation of the gradient of U
+!! in the surface layer.
+!!
+!! In section 6, the same steps are repeated but for the y direction
+!! and in section 7, a diagnostic computation of the W variance is
+!! performed.
+!!
+!! In section 8, the turbulent fluxes for the scalar variables are
+!! computed by the same way as the conservative thermodynamical variables
+!!
+!!
+!! EXTERNAL
+!! --------
+!! GX_U_M, GY_V_M, GZ_W_M : cartesian gradient operators
+!! GX_U_UW,GY_V_VW (X,Y,Z) represent the direction of the gradient
+!! _(M,U,...)_ represent the localization of the
+!! field to be derivated
+!! _(M,UW,...) represent the localization of the
+!! field derivated
+!!
+!!
+!! MXM,MXF,MYM,MYF,MZM,MZF
+!! : Shuman functions (mean operators)
+!! DXF,DYF,DZF,DZM
+!! : Shuman functions (difference operators)
+!!
+!! SUBROUTINE TRIDIAG : to compute the split implicit evolution
+!! of a variable located at a mass point
+!!
+!! SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
+!! of a variable located at a wind point
+!!
+!! FUNCTIONs ETHETA and EMOIST :
+!! allows to compute:
+!! - the coefficients for the turbulent correlation between
+!! any variable and the virtual potential temperature, of its
+!! correlations with the conservative potential temperature and
+!! the humidity conservative variable:
+!! ------- ------- -------
+!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
+!!
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : contains physical constants
+!!
+!! XG : gravity constant
+!!
+!! Module MODD_CTURB: contains the set of constants for
+!! the turbulence scheme
+!!
+!! XCMFS,XCMFB : cts for the momentum flux
+!! XCSHF : ct for the sensible heat flux
+!! XCHF : ct for the moisture flux
+!! XCTV,XCHV : cts for the T and moisture variances
+!!
+!! Module MODD_PARAMETERS
+!!
+!! JPVEXT_TURB : number of vertical external points
+!! JPHEXT : number of horizontal external points
+!!
+!!
+!! REFERENCE
+!! ---------
+!! Book 1 of documentation (Chapter: Turbulence)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original August 19, 1994
+!! Modifications: February 14, 1995 (J.Cuxart and J.Stein)
+!! Doctorization and Optimization
+!! Modifications: March 21, 1995 (J.M. Carriere)
+!! Introduction of cloud water
+!! Modifications: June 14, 1995 (J.Cuxart and J. Stein)
+!! Phi3 and Psi3 at w-point + bug in the all
+!! or nothing condens.
+!! Modifications: Sept 15, 1995 (J.Cuxart and J. Stein)
+!! Change the DP computation at the ground
+!! Modifications: October 10, 1995 (J.Cuxart and J. Stein)
+!! Psi for scal var and LES tools
+!! Modifications: November 10, 1995 (J. Stein)
+!! change the surface relations
+!! Modifications: February 20, 1995 (J. Stein) optimization
+!! Modifications: May 21, 1996 (J. Stein)
+!! bug in the vertical flux of the V wind
+!! component for explicit computation
+!! Modifications: May 21, 1996 (N. wood)
+!! modify the computation of the vertical
+!! part or the surface tangential flux
+!! Modifications: May 21, 1996 (P. Jabouille)
+!! same modification in the Y direction
+!!
+!! Modifications: Sept 17, 1996 (J. Stein) change the moist case by using
+!! Pi instead of Piref + use Atheta and Amoist
+!!
+!! Modifications: Nov 24, 1997 (V. Masson) removes the DO loops
+!! Modifications: Mar 31, 1998 (V. Masson) splits the routine TURB_VER_THERMO_FLUX
+!! Modifications: Oct 18, 2000 (V. Masson) LES computations
+!! Modifications: Dec 01, 2000 (V. Masson) conservation of energy from
+!! surface flux in 1DIM case
+!! when slopes are present
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
+!! change of YCOMMENT
+!! 2012-02 (Y. Seity) add possibility to run with reversed
+!! vertical levels
+!! Modifications July 2015 (Wim de Rooy) LHARAT switch
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!!--------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+USE MODD_CST
+USE MODD_CTURB
+USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_PARAMETERS
+USE MODD_CONF
+USE MODD_LES
+!
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_GRADIENT_M
+USE MODI_SHUMAN , ONLY : DZM, MZM, MZF
+USE MODI_TRIDIAG
+USE MODI_LES_MEAN_SUBGRID
+USE MODI_TRIDIAG_THERMO
+!
+USE MODE_IO_FIELD_WRITE, only: IO_Field_write
+USE MODE_PRANDTL
+!
+USE MODI_SECOND_MNH
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of arguments
+!
+!
+!
+INTEGER, INTENT(IN) :: KKA !near ground array index
+INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
+INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
+INTEGER, INTENT(IN) :: KRR ! number of moist var.
+INTEGER, INTENT(IN) :: KRRL ! number of liquid water var.
+INTEGER, INTENT(IN) :: KRRI ! number of ice water var.
+LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
+ ! turbulent fluxes in the syncronous FM-file
+CHARACTER(len=4), INTENT(IN) :: HTURBDIM ! dimensionality of the
+ ! turbulence scheme
+CHARACTER(len=4), INTENT(IN) :: HTOM ! type of Third Order Moment
+REAL, INTENT(IN) :: PIMPL, PEXPL ! Coef. for temporal disc.
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ, PDXX, PDYY, PDZX, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSZW ! Director Cosinus of the
+ ! normal to the ground surface
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry density * grid volum
+REAL, DIMENSION(:,:,:), INTENT(IN) :: MFMOIST ! moist mass flux dual scheme
+
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! ref. state Virtual
+ ! Potential Temperature
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHM,PSFRM ! surface fluxes at time
+! ! t - deltat
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHP,PSFRP ! surface fluxes at time
+! ! t + deltat
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM
+! Vertical wind
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM
+! potential temperature at t-Delta t
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! Mixing ratios
+ ! at t-Delta t
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! Mixing ratios
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at time t
+! In case LHARATU=TRUE, PLM already includes all stability corrections
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLEPS ! dissipative length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLOCPEXNM ! Lv(T)/Cp/Exnref at time t-1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PATHETA ! coefficients between
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PAMOIST ! s and Thetal and Rnp
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRCM ! normalized
+! 2nd-order flux s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PBETA ! buoyancy coefficient
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSQRT_TKE ! sqrt(e)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDTH_DZ ! d(th)/dz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDR_DZ ! d(rt)/dz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRED2TH3 ! 3D Redeslperger number R*2_th
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRED2R3 ! 3D Redeslperger number R*2_r
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRED2THR3 ! 3D Redeslperger number R*2_thr
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PBLL_O_E ! beta * Lk * Leps / tke
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PETHETA ! Coefficient for theta in theta_v computation
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PEMOIST ! Coefficient for r in theta_v computation
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PREDTH1 ! 1D Redelsperger number for Th
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PREDR1 ! 1D Redelsperger number for r
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPHI3 ! Prandtl number for temperature
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPSI3 ! Prandtl number for vapor
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PD ! Denominator in Prandtl numbers
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFWTH ! d(w'2th' )/dz (at flux point)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFWR ! d(w'2r' )/dz (at flux point)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFTH2 ! d(w'th'2 )/dz (at mass point)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFR2 ! d(w'r'2 )/dz (at mass point)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFTHR ! d(w'th'r')/dz (at mass point)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLP ! guess of thl at t+ deltat
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRP ! guess of r at t+ deltat
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSIGS ! Vert. part of Sigma_s at t
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+!
+REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) :: &
+ ZA, & ! work variable for wrc
+ ZFLXZ, & ! vertical flux of the treated variable
+ ZSOURCE, & ! source of evolution for the treated variable
+ ZKEFF, & ! effectif diffusion coeff = LT * SQRT( TKE )
+ ZF, & ! Flux in dTh/dt =-dF/dz (evaluated at t-1)(or rt instead of Th)
+ ZDFDDTDZ, & ! dF/d(dTh/dz)
+ ZDFDDRDZ, & ! dF/d(dr/dz)
+ Z3RDMOMENT, & ! 3 order term in flux or variance equation
+! Estimate of full level length and dissipation length scale in case LHARATU
+ PLMF, & ! estimate full level length scale from half levels (sub optimal)
+ PLEPSF ! estimate full level diss length scale from half levels (sub optimal)
+
+INTEGER :: IRESP ! Return code of FM routines
+INTEGER :: IGRID ! C-grid indicator in LFIFM file
+INTEGER :: ILENCH ! Length of comment string in LFIFM file
+INTEGER :: IKB,IKE ! I index values for the Beginning and End
+INTEGER :: IKU ! array sizes
+
+ ! mass points of the domain in the 3 direct.
+INTEGER :: I1,I2 ! For ZCOEFF allocation
+CHARACTER (LEN=100) :: YCOMMENT ! comment string in LFIFM file
+CHARACTER (LEN=16) :: YRECFM ! Name of the desired field in LFIFM file
+REAL, DIMENSION(:,:,:),ALLOCATABLE :: ZCOEFF
+ ! coefficients for the uncentred gradient
+ ! computation near the ground
+!
+REAL :: ZTIME1, ZTIME2
+!
+LOGICAL :: GUSERV ! flag to use water
+LOGICAL :: GFTH2 ! flag to use w'th'2
+LOGICAL :: GFWTH ! flag to use w'2th'
+LOGICAL :: GFR2 ! flag to use w'r'2
+LOGICAL :: GFWR ! flag to use w'2r'
+LOGICAL :: GFTHR ! flag to use w'th'r'
+TYPE(TFIELDDATA) :: TZFIELD
+!----------------------------------------------------------------------------
+!
+!* 1. PRELIMINARIES
+! -------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_CORR',0,ZHOOK_HANDLE)
+
+IKU=SIZE(PTKEM,3)
+IKB=KKA+JPVEXT_TURB*KKL
+IKE=KKU-JPVEXT_TURB*KKL
+I1=MIN(KKA+JPVEXT_TURB*KKL,KKA+JPVEXT_TURB*KKL+2*KKL)
+I2=MAX(KKA+JPVEXT_TURB*KKL,KKA+JPVEXT_TURB*KKL+2*KKL)
+
+ALLOCATE(ZCOEFF(SIZE(PDZZ,1),SIZE(PDZZ,2),I1:I2))
+!
+GUSERV = (KRR/=0)
+!
+! compute the coefficients for the uncentred gradient computation near the
+! ground
+ZCOEFF(:,:,IKB+2*KKL)= - PDZZ(:,:,IKB+KKL) / &
+ ( (PDZZ(:,:,IKB+2*KKL)+PDZZ(:,:,IKB+KKL)) * PDZZ(:,:,IKB+2*KKL) )
+ZCOEFF(:,:,IKB+KKL)= (PDZZ(:,:,IKB+2*KKL)+PDZZ(:,:,IKB+KKL)) / &
+ ( PDZZ(:,:,IKB+KKL) * PDZZ(:,:,IKB+2*KKL) )
+ZCOEFF(:,:,IKB)= - (PDZZ(:,:,IKB+2*KKL)+2.*PDZZ(:,:,IKB+KKL)) / &
+ ( (PDZZ(:,:,IKB+2*KKL)+PDZZ(:,:,IKB+KKL)) * PDZZ(:,:,IKB+KKL) )
+!
+!
+IF (LHARAT) THEN
+PLMF=MZF(PLM, KKA, KKU, KKL)
+PLEPSF=PLMF
+! function MZF produces -999 for level IKU (82 for 80 levels)
+! so put these to normal value as this level (82) is indeed calculated
+PLMF(:,:,IKU)=0.001
+PLEPSF(:,:,IKU)=0.001
+ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
+ELSE
+ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
+ENDIF
+!
+
+!
+! Flags for 3rd order quantities
+!
+GFTH2 = .FALSE.
+GFR2 = .FALSE.
+GFTHR = .FALSE.
+GFWTH = .FALSE.
+GFWR = .FALSE.
+!
+IF (HTOM/='NONE') THEN
+ GFTH2 = ANY(PFTH2/=0.)
+ GFR2 = ANY(PFR2 /=0.) .AND. GUSERV
+ GFTHR = ANY(PFTHR/=0.) .AND. GUSERV
+ GFWTH = ANY(PFWTH/=0.)
+ GFWR = ANY(PFWR /=0.) .AND. GUSERV
+END IF
+!----------------------------------------------------------------------------
+!
+!
+!* 4. TURBULENT CORRELATIONS : <THl THl>, <THl Rnp>, <Rnp Rnp>
+! --------------------------------------------------------
+!
+!
+!* 4.2 <THl THl>
+!
+! Compute the turbulent variance F and F' at time t-dt.
+!
+IF (LHARAT) THEN
+ ZF (:,:,:) = PLMF*PLEPSF*MZF(PDTH_DZ**2, KKA, KKU, KKL)
+ELSE
+ ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(PPHI3*PDTH_DZ**2, KKA, KKU, KKL)
+ENDIF
+ ZDFDDTDZ(:,:,:) = 0. ! this term, because of discretization, is treated separately
+ !
+ ! Effect of 3rd order terms in temperature flux (at mass point)
+ !
+ ! d(w'th'2)/dz
+ IF (GFTH2) THEN
+ ZF = ZF + M3_TH2_WTH2(KKA,KKU,KKL,PREDTH1,PREDR1,PD,PLEPS,&
+ & PSQRT_TKE) * PFTH2
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_WTH2_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,&
+ & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA) * PFTH2
+ END IF
+ !
+ ! d(w'2th')/dz
+ IF (GFWTH) THEN
+ ZF = ZF + M3_TH2_W2TH(KKA,KKU,KKL,PREDTH1,PREDR1,PD,PDTH_DZ,&
+ & PLM,PLEPS,PTKEM) * MZF(PFWTH, KKA, KKU, KKL)
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_W2TH_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,PD,&
+ & PLM,PLEPS,PTKEM,GUSERV) * MZF(PFWTH, KKA, KKU, KKL)
+ END IF
+ !
+ IF (KRR/=0) THEN
+ ! d(w'r'2)/dz
+ IF (GFR2) THEN
+ ZF = ZF + M3_TH2_WR2(KKA,KKU,KKL,PD,PLEPS,PSQRT_TKE,PBLL_O_E,&
+ & PEMOIST,PDTH_DZ) * PFR2
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_WR2_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,PD,&
+ & PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST,PDTH_DZ) * PFR2
+ END IF
+ !
+ ! d(w'2r')/dz
+ IF (GFWR) THEN
+ ZF = ZF + M3_TH2_W2R(KKA,KKU,KKL,PD,PLM,PLEPS,PTKEM,PBLL_O_E,&
+ & PEMOIST,PDTH_DZ) * MZF(PFWR, KKA, KKU, KKL)
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_W2R_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,PD,&
+ & PLM,PLEPS,PTKEM,PBLL_O_E,PEMOIST,PDTH_DZ) * MZF(PFWR, KKA, KKU, KKL)
+ END IF
+ !
+ ! d(w'th'r')/dz
+ IF (GFTHR) THEN
+ ZF = ZF + M3_TH2_WTHR(KKA,KKU,KKL,PREDR1,PD,PLEPS,PSQRT_TKE,&
+ & PBLL_O_E,PEMOIST,PDTH_DZ) * PFTHR
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_TH2_WTHR_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,&
+ & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST,PDTH_DZ) * PFTHR
+ END IF
+
+ END IF
+ !
+ ZFLXZ(:,:,:) = ZF &
+ ! + PIMPL * XCTV*PLM*PLEPS &
+ ! *MZF(D_PHI3DTDZ2_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,PDTH_DZ,HTURBDIM,GUSERV) &
+ ! *DZM(PTHLP - PTHLM, KKA, KKU, KKL) / PDZZ ) &
+ + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM, KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
+ !
+ ! special case near the ground ( uncentred gradient )
+ IF (LHARAT) THEN
+ ZFLXZ(:,:,IKB) = PLMF(:,:,IKB) &
+ * PLEPSF(:,:,IKB) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB ) )**2 &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB ) )**2 &
+ )
+ ELSE
+ ZFLXZ(:,:,IKB) = XCTV * PPHI3(:,:,IKB+KKL) * PLM(:,:,IKB) &
+ * PLEPS(:,:,IKB) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB ) )**2 &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB ) )**2 &
+ )
+ ENDIF
+ !
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ !
+ ZFLXZ = MAX(0., ZFLXZ)
+ !
+ IF (KRRL > 0) THEN
+ PSIGS(:,:,:) = ZFLXZ(:,:,:) * PATHETA(:,:,:)**2
+ END IF
+ !
+ !
+ ! stores <THl THl>
+ IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+ TZFIELD%CMNHNAME = 'THL_VVAR'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'THL_VVAR'
+ TZFIELD%CUNITS = 'K2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_THL_VVAR'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
+ END IF
+!
+! and we store in LES configuration
+!
+ IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID(ZFLXZ, X_LES_SUBGRID_Thl2 )
+ CALL LES_MEAN_SUBGRID(MZF(PWM, KKA, KKU, KKL)*ZFLXZ, X_LES_RES_W_SBG_Thl2 )
+ CALL LES_MEAN_SUBGRID(-2.*XCTD*PSQRT_TKE*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Thl2 )
+ CALL LES_MEAN_SUBGRID(PETHETA*ZFLXZ, X_LES_SUBGRID_ThlThv )
+ CALL LES_MEAN_SUBGRID(-XA3*PBETA*PETHETA*ZFLXZ, X_LES_SUBGRID_ThlPz, .TRUE. )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+!
+ IF ( KRR /= 0 ) THEN
+!
+!* 4.3 <THl Rnp>
+!
+!
+ ! Compute the turbulent variance F and F' at time t-dt.
+IF (LHARAT) THEN
+ ZF (:,:,:) = PLMF*PLEPSF*MZF(PDTH_DZ*PDR_DZ, KKA, KKU, KKL)
+ELSE
+ ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(0.5*(PPHI3+PPSI3)*PDTH_DZ*PDR_DZ, KKA, KKU, KKL)
+ENDIF
+ ZDFDDTDZ(:,:,:) = 0. ! this term, because of discretization, is treated separately
+ ZDFDDRDZ(:,:,:) = 0. ! this term, because of discretization, is treated separately
+ !
+ ! Effect of 3rd order terms in temperature flux (at mass point)
+ !
+ ! d(w'th'2)/dz
+ IF (GFTH2) THEN
+ ZF = ZF + M3_THR_WTH2(KKA,KKU,KKL,PREDR1,PD,PLEPS,PSQRT_TKE,&
+ & PBLL_O_E,PETHETA,PDR_DZ) * PFTH2
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_WTH2_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,&
+ & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA,PDR_DZ) * PFTH2
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_WTH2_O_DDRDZ(KKA,KKU,KKL,PREDTH1,PREDR1,&
+ & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA) * PFTH2
+ END IF
+ !
+ ! d(w'2th')/dz
+ IF (GFWTH) THEN
+ ZF = ZF + M3_THR_W2TH(KKA,KKU,KKL,PREDR1,PD,PLM,PLEPS,PTKEM,&
+ & PDR_DZ) * MZF(PFWTH, KKA, KKU, KKL)
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_W2TH_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,&
+ & PD,PLM,PLEPS,PTKEM,PBLL_O_E,PDR_DZ,PETHETA) * MZF(PFWTH, KKA, KKU, KKL)
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_W2TH_O_DDRDZ(KKA,KKU,KKL,PREDTH1,PREDR1,&
+ & PD,PLM,PLEPS,PTKEM) * MZF(PFWTH, KKA, KKU, KKL)
+ END IF
+ !
+ ! d(w'r'2)/dz
+ IF (GFR2) THEN
+ ZF = ZF + M3_THR_WR2(KKA,KKU,KKL,PREDTH1,PD,PLEPS,PSQRT_TKE,&
+ & PBLL_O_E,PEMOIST,PDTH_DZ) * PFR2
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_WR2_O_DDTDZ(KKA,KKU,KKL,PREDR1,PREDTH1,PD,&
+ & PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST) * PFR2
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_WR2_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,PD,&
+ & PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST,PDTH_DZ) * PFR2
+ END IF
+ !
+ ! d(w'2r')/dz
+ IF (GFWR) THEN
+ ZF = ZF + M3_THR_W2R(KKA,KKU,KKL,PREDTH1,PD,PLM,PLEPS,PTKEM,&
+ & PDTH_DZ) * MZF(PFWR, KKA, KKU, KKL)
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_W2R_O_DDTDZ(KKA,KKU,KKL,PREDR1,PREDTH1,PD,&
+ & PLM,PLEPS,PTKEM) * MZF(PFWR, KKA, KKU, KKL)
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_W2R_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,PD,&
+ & PLM,PLEPS,PTKEM,PBLL_O_E,PDTH_DZ,PEMOIST) * MZF(PFWR, KKA, KKU, KKL)
+ END IF
+ !
+ ! d(w'th'r')/dz
+ IF (GFTHR) THEN
+ ZF = ZF + M3_THR_WTHR(KKA,KKU,KKL,PREDTH1,PREDR1,PD,PLEPS,&
+ & PSQRT_TKE) * PFTHR
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_THR_WTHR_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,&
+ & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA) * PFTHR
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_THR_WTHR_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,&
+ & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST) * PFTHR
+ END IF
+ !
+ IF (LHARAT) THEN
+ ZFLXZ(:,:,:) = ZF &
+ + PIMPL * PLMF*PLEPSF*0.5 &
+ * MZF(( 2. &
+ ) *PDR_DZ *DZM(PTHLP - PTHLM, KKA, KKU, KKL ) / PDZZ &
+ +( 2. &
+ ) *PDTH_DZ *DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ &
+ , KKA, KKU, KKL) &
+ + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM(:,:,:), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
+ + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
+ ELSE
+ ZFLXZ(:,:,:) = ZF &
+ + PIMPL * XCTV*PLM*PLEPS*0.5 &
+ * MZF(( D_PHI3DTDZ_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV) & ! d(phi3*dthdz)/ddthdz term
+ +D_PSI3DTDZ_O_DDTDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV) & ! d(psi3*dthdz)/ddthdz term
+ ) *PDR_DZ *DZM(PTHLP - PTHLM, KKA, KKU, KKL) / PDZZ &
+ +( D_PHI3DRDZ_O_DDRDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV) & ! d(phi3*drdz )/ddrdz term
+ +D_PSI3DRDZ_O_DDRDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV) & ! d(psi3*drdz )/ddrdz term
+ ) *PDTH_DZ *DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ &
+ , KKA, KKU, KKL) &
+ + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM(:,:,:), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
+ + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
+ ENDIF
+ !
+ ! special case near the ground ( uncentred gradient )
+ IF (LHARAT) THEN
+ ZFLXZ(:,:,IKB) = &
+ (1. ) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB )) &
+ *( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
+ +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
+ +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 )) &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB )) &
+ *( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL ) &
+ +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PRP(:,:,IKB )) &
+ )
+ ELSE
+ ZFLXZ(:,:,IKB) = &
+ (XCHT1 * PPHI3(:,:,IKB+KKL) + XCHT2 * PPSI3(:,:,IKB+KKL)) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB )) &
+ *( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
+ +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
+ +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 )) &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB )) &
+ *( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL ) &
+ +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PRP(:,:,IKB )) &
+ )
+ ENDIF
+ !
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ !
+ IF ( KRRL > 0 ) THEN
+!
+!
+! NB PATHETA is -b in Chaboureau Bechtold 2002 which explains the + sign here
+
+ PSIGS(:,:,:) = PSIGS(:,:,:) + &
+ 2. * PATHETA(:,:,:) * PAMOIST(:,:,:) * ZFLXZ(:,:,:)
+ END IF
+ ! stores <THl Rnp>
+ IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+ TZFIELD%CMNHNAME = 'THLRCONS_VCOR'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'THLRCONS_VCOR'
+ TZFIELD%CUNITS = 'K kg kg-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_THLRCONS_VCOR'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
+ END IF
+!
+! and we store in LES configuration
+!
+IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID(ZFLXZ, X_LES_SUBGRID_THlRt )
+ CALL LES_MEAN_SUBGRID(MZF(PWM, KKA, KKU, KKL)*ZFLXZ, X_LES_RES_W_SBG_ThlRt )
+ CALL LES_MEAN_SUBGRID(-2.*XCTD*PSQRT_TKE*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_ThlRt )
+ CALL LES_MEAN_SUBGRID(PETHETA*ZFLXZ, X_LES_SUBGRID_RtThv )
+ CALL LES_MEAN_SUBGRID(-XA3*PBETA*PETHETA*ZFLXZ, X_LES_SUBGRID_RtPz, .TRUE. )
+ CALL LES_MEAN_SUBGRID(PEMOIST*ZFLXZ, X_LES_SUBGRID_ThlThv , .TRUE. )
+ CALL LES_MEAN_SUBGRID(-XA3*PBETA*PEMOIST*ZFLXZ, X_LES_SUBGRID_ThlPz, .TRUE. )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+!
+!
+!* 4.4 <Rnp Rnp>
+!
+!
+ ! Compute the turbulent variance F and F' at time t-dt.
+IF (LHARAT) THEN
+ ZF (:,:,:) = PLMF*PLEPSF*MZF(PDR_DZ**2, KKA, KKU, KKL)
+ ELSE
+ ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(PPSI3*PDR_DZ**2, KKA, KKU, KKL)
+ENDIF
+ ZDFDDRDZ(:,:,:) = 0. ! this term, because of discretization, is treated separately
+ !
+ ! Effect of 3rd order terms in temperature flux (at mass point)
+ !
+ ! d(w'r'2)/dz
+ IF (GFR2) THEN
+ ZF = ZF + M3_R2_WR2(KKA,KKU,KKL,PREDR1,PREDTH1,PD,PLEPS,&
+ & PSQRT_TKE) * PFR2
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_WR2_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,&
+ & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PEMOIST) * PFR2
+ END IF
+ !
+ ! d(w'2r')/dz
+ IF (GFWR) THEN
+ ZF = ZF + M3_R2_W2R(KKA,KKU,KKL,PREDR1,PREDTH1,PD,PDR_DZ,&
+ & PLM,PLEPS,PTKEM) * MZF(PFWR, KKA, KKU, KKL)
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_W2R_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,&
+ & PD,PLM,PLEPS,PTKEM,GUSERV) * MZF(PFWR, KKA, KKU, KKL)
+ END IF
+ !
+ IF (KRR/=0) THEN
+ ! d(w'r'2)/dz
+ IF (GFTH2) THEN
+ ZF = ZF + M3_R2_WTH2(KKA,KKU,KKL,PD,PLEPS,PSQRT_TKE,&
+ & PBLL_O_E,PETHETA,PDR_DZ) * PFTH2
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_WTH2_O_DDRDZ(KKA,KKU,KKL,PREDR1,&
+ & PREDTH1,PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA,PDR_DZ) * PFTH2
+ END IF
+ !
+ ! d(w'2r')/dz
+ IF (GFWTH) THEN
+ ZF = ZF + M3_R2_W2TH(KKA,KKU,KKL,PD,PLM,PLEPS,PTKEM,&
+ & PBLL_O_E,PETHETA,PDR_DZ) * MZF(PFWTH, KKA, KKU, KKL)
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_W2TH_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,&
+ & PD,PLM,PLEPS,PTKEM,PBLL_O_E,PETHETA,PDR_DZ) * MZF(PFWTH, KKA, KKU, KKL)
+ END IF
+ !
+ ! d(w'th'r')/dz
+ IF (GFTHR) THEN
+ ZF = ZF + M3_R2_WTHR(KKA,KKU,KKL,PREDTH1,PD,PLEPS,&
+ & PSQRT_TKE,PBLL_O_E,PETHETA,PDR_DZ) * PFTHR
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_WTHR_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,&
+ & PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA,PDR_DZ) * PFTHR
+ END IF
+
+ END IF
+ !
+ IF (LHARAT) THEN
+ ZFLXZ(:,:,:) = ZF &
+ + PIMPL * PLMF*PLEPSF &
+ *MZF(2.*PDR_DZ* &
+ DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
+ + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
+ ELSE
+ ZFLXZ(:,:,:) = ZF &
+ + PIMPL * XCTV*PLM*PLEPS &
+ *MZF(D_PSI3DRDZ2_O_DDRDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,PDR_DZ,HTURBDIM,GUSERV) &
+ *DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
+ + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
+ ENDIF
+ !
+ ! special case near the ground ( uncentred gradient )
+ IF (LHARAT) THEN
+ ZFLXZ(:,:,IKB) = PLMF(:,:,IKB) &
+ * PLEPSF(:,:,IKB) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
+ +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
+ +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 ))**2 &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PRP(:,:,IKB ))**2 &
+ )
+ ELSE
+ ZFLXZ(:,:,IKB) = XCHV * PPSI3(:,:,IKB+KKL) * PLM(:,:,IKB) &
+ * PLEPS(:,:,IKB) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
+ +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
+ +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 ))**2 &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PRP(:,:,IKB ))**2 &
+ )
+ ENDIF
+ !
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ !
+ IF ( KRRL > 0 ) THEN
+ PSIGS(:,:,:) = PSIGS(:,:,:) + PAMOIST(:,:,:) **2 * ZFLXZ(:,:,:)
+ END IF
+ ! stores <Rnp Rnp>
+ IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+ TZFIELD%CMNHNAME = 'RTOT_VVAR'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'RTOT_VVAR'
+ TZFIELD%CUNITS = 'kg2 kg-2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_RTOT_VVAR'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
+ END IF
+ !
+ ! and we store in LES configuration
+ !
+ IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID(ZFLXZ, X_LES_SUBGRID_Rt2 )
+ CALL LES_MEAN_SUBGRID(MZF(PWM, KKA, KKU, KKL)*ZFLXZ, X_LES_RES_W_SBG_Rt2 )
+ CALL LES_MEAN_SUBGRID(PEMOIST*ZFLXZ, X_LES_SUBGRID_RtThv , .TRUE. )
+ CALL LES_MEAN_SUBGRID(-XA3*PBETA*PEMOIST*ZFLXZ, X_LES_SUBGRID_RtPz, .TRUE. )
+ CALL LES_MEAN_SUBGRID(-2.*XCTD*PSQRT_TKE*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Rt2 )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+ !
+ END IF ! end if KRR ne 0
+!
+!
+! 4.5 Vertical part of Sigma_s
+!
+ IF ( KRRL > 0 ) THEN
+ ! Extrapolate PSIGS at the ground and at the top
+ PSIGS(:,:,KKA) = PSIGS(:,:,IKB)
+ PSIGS(:,:,KKU) = PSIGS(:,:,IKE)
+ PSIGS(:,:,:) = MAX (PSIGS(:,:,:) , 0.)
+ PSIGS(:,:,:) = SQRT(PSIGS(:,:,:))
+ END IF
+
+!
+! 4.6 Deallocate
+!
+ DEALLOCATE(ZCOEFF)
+!----------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_CORR',1,ZHOOK_HANDLE)
+END SUBROUTINE TURB_VER_THERMO_CORR
+END MODULE MODE_TURB_VER_THERMO_CORR
diff --git a/src/mesonh/turb/mode_turb_ver_thermo_flux.f90 b/src/mesonh/turb/mode_turb_ver_thermo_flux.f90
new file mode 100644
index 000000000..7a79162a4
--- /dev/null
+++ b/src/mesonh/turb/mode_turb_ver_thermo_flux.f90
@@ -0,0 +1,833 @@
+!MNH_LIC Copyright 1994-2021 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 MODE_TURB_VER_THERMO_FLUX
+IMPLICIT NONE
+CONTAINS
+
+SUBROUTINE TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI, &
+ OTURB_FLX,HTURBDIM,HTOM, &
+ PIMPL,PEXPL, &
+ PTSTEP, &
+ TPFILE, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ, &
+ PRHODJ,PTHVREF, &
+ PSFTHM,PSFRM,PSFTHP,PSFRP, &
+ PWM,PTHLM,PRM,PSVM, &
+ PTKEM,PLM,PLEPS, &
+ PLOCPEXNM,PATHETA,PAMOIST,PSRCM,PFRAC_ICE, &
+ PBETA, PSQRT_TKE, PDTH_DZ, PDR_DZ, PRED2TH3, &
+ PRED2R3, PRED2THR3, PBLL_O_E, PETHETA, &
+ PEMOIST, PREDTH1, PREDR1, PPHI3, PPSI3, PD, &
+ PFWTH,PFWR,PFTH2,PFR2,PFTHR,MFMOIST,PBL_DEPTH,&
+ PWTHV,PRTHLS,PRRS,PTHLP,PRP,PTP,PWTH,PWRC )
+! ###############################################################
+!
+!
+!!**** *TURB_VER_THERMO_FLUX* -compute the source terms due to the vertical turbulent
+!! fluxes.
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute the vertical turbulent
+! fluxes of the evolutive variables and give back the source
+! terms to the main program. In the case of large horizontal meshes,
+! the divergence of these vertical turbulent fluxes represent the whole
+! effect of the turbulence but when the three-dimensionnal version of
+! the turbulence scheme is activated (CTURBDIM="3DIM"), these divergences
+! are completed in the next routine TURB_HOR.
+! An arbitrary degree of implicitness has been implemented for the
+! temporal treatment of these diffusion terms.
+! The vertical boundary conditions are as follows:
+! * at the bottom, the surface fluxes are prescribed at the same
+! as the other turbulent fluxes
+! * at the top, the turbulent fluxes are set to 0.
+! It should be noted that the condensation has been implicitely included
+! in this turbulence scheme by using conservative variables and computing
+! the subgrid variance of a statistical variable s indicating the presence
+! or not of condensation in a given mesh.
+!
+!!** METHOD
+!! ------
+!! 1D type calculations are made;
+!! The vertical turbulent fluxes are computed in an off-centered
+!! implicit scheme (a Crank-Nicholson type with coefficients different
+!! than 0.5), which allows to vary the degree of implicitness of the
+!! formulation.
+!! The different prognostic variables are treated one by one.
+!! The contributions of each turbulent fluxes are cumulated into the
+!! tendency PRvarS, and into the dynamic and thermal production of
+!! TKE if necessary.
+!!
+!! In section 2 and 3, the thermodynamical fields are considered.
+!! Only the turbulent fluxes of the conservative variables
+!! (Thetal and Rnp stored in PRx(:,:,:,1)) are computed.
+!! Note that the turbulent fluxes at the vertical
+!! boundaries are given either by the soil scheme for the surface one
+!! ( at the same instant as the others fluxes) and equal to 0 at the
+!! top of the model. The thermal production is computed by vertically
+!! averaging the turbulent flux and multiply this flux at the mass point by
+!! a function ETHETA or EMOIST, which preform the transformation from the
+!! conservative variables to the virtual potential temperature.
+!!
+!! In section 4, the variance of the statistical variable
+!! s indicating presence or not of condensation, is determined in function
+!! of the turbulent moments of the conservative variables and its
+!! squarred root is stored in PSIGS. This information will be completed in
+!! the horizontal turbulence if the turbulence dimensionality is not
+!! equal to "1DIM".
+!!
+!! In section 5, the x component of the stress tensor is computed.
+!! The surface flux <u'w'> is computed from the value of the surface
+!! fluxes computed in axes linked to the orography ( i", j" , k"):
+!! i" is parallel to the surface and in the direction of the maximum
+!! slope
+!! j" is also parallel to the surface and in the normal direction of
+!! the maximum slope
+!! k" is the normal to the surface
+!! In order to prevent numerical instability, the implicit scheme has
+!! been extended to the surface flux regarding to its dependence in
+!! function of U. The dependence in function of the other components
+!! introduced by the different rotations is only explicit.
+!! The turbulent fluxes are used to compute the dynamic production of
+!! TKE. For the last TKE level ( located at PDZZ(:,:,IKB)/2 from the
+!! ground), an harmonic extrapolation from the dynamic production at
+!! PDZZ(:,:,IKB) is used to avoid an evaluation of the gradient of U
+!! in the surface layer.
+!!
+!! In section 6, the same steps are repeated but for the y direction
+!! and in section 7, a diagnostic computation of the W variance is
+!! performed.
+!!
+!! In section 8, the turbulent fluxes for the scalar variables are
+!! computed by the same way as the conservative thermodynamical variables
+!!
+!!
+!! EXTERNAL
+!! --------
+!! GX_U_M, GY_V_M, GZ_W_M : cartesian gradient operators
+!! GX_U_UW,GY_V_VW (X,Y,Z) represent the direction of the gradient
+!! _(M,U,...)_ represent the localization of the
+!! field to be derivated
+!! _(M,UW,...) represent the localization of the
+!! field derivated
+!!
+!!
+!! MXM,MXF,MYM,MYF,MZM,MZF
+!! : Shuman functions (mean operators)
+!! DXF,DYF,DZF,DZM
+!! : Shuman functions (difference operators)
+!!
+!! SUBROUTINE TRIDIAG : to compute the split implicit evolution
+!! of a variable located at a mass point
+!!
+!! SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
+!! of a variable located at a wind point
+!!
+!! FUNCTIONs ETHETA and EMOIST :
+!! allows to compute:
+!! - the coefficients for the turbulent correlation between
+!! any variable and the virtual potential temperature, of its
+!! correlations with the conservative potential temperature and
+!! the humidity conservative variable:
+!! ------- ------- -------
+!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
+!!
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : contains physical constants
+!!
+!! XG : gravity constant
+!!
+!! Module MODD_CTURB: contains the set of constants for
+!! the turbulence scheme
+!!
+!! XCMFS,XCMFB : cts for the momentum flux
+!! XCSHF : ct for the sensible heat flux
+!! XCHF : ct for the moisture flux
+!! XCTV,XCHV : cts for the T and moisture variances
+!!
+!! Module MODD_PARAMETERS
+!!
+!! JPVEXT_TURB : number of vertical external points
+!! JPHEXT : number of horizontal external points
+!!
+!!
+!! REFERENCE
+!! ---------
+!! Book 1 of documentation (Chapter: Turbulence)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original August 19, 1994
+!! Modifications: February 14, 1995 (J.Cuxart and J.Stein)
+!! Doctorization and Optimization
+!! Modifications: March 21, 1995 (J.M. Carriere)
+!! Introduction of cloud water
+!! Modifications: June 14, 1995 (J.Cuxart and J. Stein)
+!! Phi3 and Psi3 at w-point + bug in the all
+!! or nothing condens.
+!! Modifications: Sept 15, 1995 (J.Cuxart and J. Stein)
+!! Change the DP computation at the ground
+!! Modifications: October 10, 1995 (J.Cuxart and J. Stein)
+!! Psi for scal var and LES tools
+!! Modifications: November 10, 1995 (J. Stein)
+!! change the surface relations
+!! Modifications: February 20, 1995 (J. Stein) optimization
+!! Modifications: May 21, 1996 (J. Stein)
+!! bug in the vertical flux of the V wind
+!! component for explicit computation
+!! Modifications: May 21, 1996 (N. wood)
+!! modify the computation of the vertical
+!! part or the surface tangential flux
+!! Modifications: May 21, 1996 (P. Jabouille)
+!! same modification in the Y direction
+!!
+!! Modifications: Sept 17, 1996 (J. Stein) change the moist case by using
+!! Pi instead of Piref + use Atheta and Amoist
+!!
+!! Modifications: Nov 24, 1997 (V. Masson) removes the DO loops
+!! Modifications: Mar 31, 1998 (V. Masson) splits the routine TURB_VER_THERMO_FLUX
+!! Modifications: Oct 18, 2000 (V. Masson) LES computations
+!! Modifications: Dec 01, 2000 (V. Masson) conservation of energy from
+!! surface flux in 1DIM case
+!! when slopes are present
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! Feb 20, 2003 (JP Pinty) Add PFRAC_ICE
+!! May 20, 2003 (JP Pinty) Correction of ETHETA
+!! and EMOIST calls
+!! July 2005 (S. Tomas, V. Masson)
+!! Add 3rd order moments
+!! and implicitation of PHI3 and PSI3
+!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
+!! change of YCOMMENT
+!! 2012-02 (Y. Seity) add possibility to run with reversed
+!! vertical levels
+!! Modifications July 2015 (Wim de Rooy) LHARAT switch
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! 2021 (D. Ricard) last version of HGRAD turbulence scheme
+!! Leronard terms instead of Reynolds terms
+!! applied to vertical fluxes of r_np and Thl
+!! for implicit version of turbulence scheme
+!! corrections and cleaning
+!! June 2020 (B. Vie) Patch preventing negative rc and ri in 2.3 and 3.3
+!! JL Redelsperger : 03/2021: Ocean and Autocoupling O-A LES Cases
+!! Sfc flux shape for LDEEPOC Case
+!!--------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
+USE MODD_CST
+USE MODD_CTURB
+USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
+USE MODD_GRID_n, ONLY: XZS, XXHAT, XYHAT
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_METRICS_n, ONLY: XDXX, XDYY, XDZX, XDZY, XDZZ
+USE MODD_PARAMETERS
+USE MODD_TURB_n, ONLY: LHGRAD, XCOEFHGRADTHL, XCOEFHGRADRM, XALTHGRAD, XCLDTHOLD
+USE MODD_CONF
+USE MODD_LES
+USE MODD_DIM_n
+USE MODD_DYN_n, ONLY: LOCEAN
+USE MODD_OCEANH
+USE MODD_REF, ONLY: LCOUPLES
+USE MODD_TURB_n
+USE MODD_FRC
+!
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_GRADIENT_M
+USE MODI_SHUMAN , ONLY : DZF, DZM, MZF, MZM
+USE MODI_TRIDIAG
+USE MODI_LES_MEAN_SUBGRID
+USE MODI_TRIDIAG_THERMO
+USE MODI_TM06_H
+!
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
+USE MODE_PRANDTL
+!
+USE MODI_SECOND_MNH
+USE MODE_ll
+USE MODE_GATHER_ll
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of arguments
+!
+!
+!
+INTEGER, INTENT(IN) :: KKA !near ground array index
+INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
+INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
+INTEGER, INTENT(IN) :: KRR ! number of moist var.
+INTEGER, INTENT(IN) :: KRRL ! number of liquid water var.
+INTEGER, INTENT(IN) :: KRRI ! number of ice water var.
+LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
+ ! turbulent fluxes in the syncronous FM-file
+CHARACTER(len=4), INTENT(IN) :: HTURBDIM ! dimensionality of the
+ ! turbulence scheme
+CHARACTER(len=4), INTENT(IN) :: HTOM ! type of Third Order Moment
+REAL, INTENT(IN) :: PIMPL, PEXPL ! Coef. for temporal disc.
+REAL, INTENT(IN) :: PTSTEP ! Double Time Step
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ, PDXX, PDYY, PDZX, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSZW ! Director Cosinus of the
+ ! normal to the ground surface
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! altitudes
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry density * grid volum
+REAL, DIMENSION(:,:,:), INTENT(IN) :: MFMOIST ! moist mass flux dual scheme
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! ref. state Virtual
+ ! Potential Temperature
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHM,PSFRM ! surface fluxes at time
+! ! t - deltat
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHP,PSFRP ! surface fluxes at time
+! ! t + deltat
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM
+! Vertical wind
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM
+! potential temperature at t-Delta t
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! Mixing ratios
+ ! at t-Delta t
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! Mixing ratios
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at time t
+!
+! In case LHARAT=TRUE, PLM already includes all stability corrections
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLEPS ! dissipative length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLOCPEXNM ! Lv(T)/Cp/Exnref at time t-1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PATHETA ! coefficients between
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PAMOIST ! s and Thetal and Rnp
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRCM ! normalized
+! 2nd-order flux s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFRAC_ICE ! ri fraction of rc+ri
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PBETA ! buoyancy coefficient
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSQRT_TKE ! sqrt(e)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDTH_DZ ! d(th)/dz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDR_DZ ! d(rt)/dz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRED2TH3 ! 3D Redeslperger number R*2_th
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRED2R3 ! 3D Redeslperger number R*2_r
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRED2THR3 ! 3D Redeslperger number R*2_thr
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PBLL_O_E ! beta * Lk * Leps / tke
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PETHETA ! Coefficient for theta in theta_v computation
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PEMOIST ! Coefficient for r in theta_v computation
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PREDTH1 ! 1D Redelsperger number for Th
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PREDR1 ! 1D Redelsperger number for r
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPHI3 ! Prandtl number for temperature
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPSI3 ! Prandtl number for vapor
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PD ! Denominator in Prandtl numbers
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFWTH ! d(w'2th' )/dz (at flux point)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFWR ! d(w'2r' )/dz (at flux point)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFTH2 ! d(w'th'2 )/dz (at mass point)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFR2 ! d(w'r'2 )/dz (at mass point)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFTHR ! d(w'th'r')/dz (at mass point)
+REAL, DIMENSION(:,:), INTENT(INOUT):: PBL_DEPTH ! BL depth
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PWTHV ! buoyancy flux
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHLS ! cumulated source for theta
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS ! cumulated source for rt
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PTHLP ! guess of thl at t+ deltat
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRP ! guess of r at t+ deltat
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PTP ! Dynamic and thermal
+ ! TKE production terms
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PWTH ! heat flux
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PWRC ! cloud water flux
+!
+!
+!* 0.2 declaration of local variables
+!
+!
+REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) :: &
+ ZA, & ! work variable for wrc or LES computation
+ ZFLXZ, & ! vertical flux of the treated variable
+ ZSOURCE, & ! source of evolution for the treated variable
+ ZKEFF, & ! effectif diffusion coeff = LT * SQRT( TKE )
+ ZF, & ! Flux in dTh/dt =-dF/dz (evaluated at t-1)(or rt instead of Th)
+ ZDFDDTDZ, & ! dF/d(dTh/dz)
+ ZDFDDRDZ, & ! dF/d(dr/dz)
+ Z3RDMOMENT ! 3 order term in flux or variance equation
+INTEGER :: IRESP ! Return code of FM routines
+INTEGER :: IGRID ! C-grid indicator in LFIFM file
+INTEGER :: ILENCH ! Length of comment string in LFIFM file
+INTEGER :: IKB,IKE ! I index values for the Beginning and End
+ ! mass points of the domain in the 3 direct.
+INTEGER :: IKT ! array size in k direction
+INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
+CHARACTER (LEN=100) :: YCOMMENT ! comment string in LFIFM file
+CHARACTER (LEN=16) :: YRECFM ! Name of the desired field in LFIFM file
+!
+REAL :: ZTIME1, ZTIME2
+!
+INTEGER :: JK
+LOGICAL :: GUSERV ! flag to use water
+LOGICAL :: GFTH2 ! flag to use w'th'2
+LOGICAL :: GFWTH ! flag to use w'2th'
+LOGICAL :: GFR2 ! flag to use w'r'2
+LOGICAL :: GFWR ! flag to use w'2r'
+LOGICAL :: GFTHR ! flag to use w'th'r'
+TYPE(TFIELDDATA) :: TZFIELD
+!----------------------------------------------------------------------------
+!
+!* 1. PRELIMINARIES
+! -------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_FLUX',0,ZHOOK_HANDLE)
+IKT =SIZE(PTHLM,3)
+IKTE =IKT-JPVEXT_TURB
+IKTB =1+JPVEXT_TURB
+IKB=KKA+JPVEXT_TURB*KKL
+IKE=KKU-JPVEXT_TURB*KKL
+!
+GUSERV = (KRR/=0)
+!
+! compute the coefficients for the uncentred gradient computation near the
+! ground
+!
+IF (LHARAT) THEN
+! LHARAT so TKE and length scales at half levels!
+ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) +50.*MFMOIST(:,:,:)
+ELSE
+ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
+ENDIF
+!
+!
+! Flags for 3rd order quantities
+!
+GFTH2 = .FALSE.
+GFR2 = .FALSE.
+GFTHR = .FALSE.
+GFWTH = .FALSE.
+GFWR = .FALSE.
+!
+IF (HTOM/='NONE') THEN
+ GFTH2 = ANY(PFTH2/=0.)
+ GFR2 = ANY(PFR2 /=0.) .AND. GUSERV
+ GFTHR = ANY(PFTHR/=0.) .AND. GUSERV
+ GFWTH = ANY(PFWTH/=0.)
+ GFWR = ANY(PFWR /=0.) .AND. GUSERV
+END IF
+!----------------------------------------------------------------------------
+!
+!* 2. SOURCES OF CONSERVATIVE POTENTIAL TEMPERATURE AND
+! PARTIAL THERMAL PRODUCTION
+! ---------------------------------------------------------------
+!
+!* 2.1 Splitted value for cons. potential temperature at t+deltat
+!
+! Compute the turbulent flux F and F' at time t-dt.
+!
+IF (LHARAT) THEN
+ZF (:,:,:) = -ZKEFF*DZM(PTHLM, KKA, KKU, KKL)/PDZZ
+ZDFDDTDZ(:,:,:) = -ZKEFF
+ELSE
+ZF (:,:,:) = -XCSHF*PPHI3*ZKEFF*DZM(PTHLM, KKA, KKU, KKL)/PDZZ
+ZDFDDTDZ(:,:,:) = -XCSHF*ZKEFF*D_PHI3DTDZ_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV)
+ENDIF
+!
+! Effect of 3rd order terms in temperature flux (at flux point)
+!
+! d(w'2th')/dz
+IF (GFWTH) THEN
+ Z3RDMOMENT= M3_WTH_W2TH(KKA,KKU,KKL,PREDTH1,PREDR1,PD,ZKEFF,PTKEM)
+!
+ ZF = ZF + Z3RDMOMENT * PFWTH
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_W2TH_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,&
+ & PD,PBLL_O_E,PETHETA,ZKEFF,PTKEM) * PFWTH
+END IF
+!
+! d(w'th'2)/dz
+IF (GFTH2) THEN
+ Z3RDMOMENT= M3_WTH_WTH2(PREDTH1,PREDR1,PD,PBLL_O_E,PETHETA)
+!
+ ZF = ZF + Z3RDMOMENT * MZM(PFTH2, KKA, KKU, KKL)
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_WTH2_O_DDTDZ(Z3RDMOMENT,PREDTH1,PREDR1,&
+ & PD,PBLL_O_E,PETHETA) * MZM(PFTH2, KKA, KKU, KKL)
+END IF
+!
+! d(w'2r')/dz
+IF (GFWR) THEN
+ ZF = ZF + M3_WTH_W2R(KKA,KKU,KKL,PREDTH1,PREDR1,PD,ZKEFF,&
+ & PTKEM,PBLL_O_E,PEMOIST,PDTH_DZ) * PFWR
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_W2R_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,&
+ & PD,ZKEFF,PTKEM,PBLL_O_E,PEMOIST) * PFWR
+END IF
+!
+! d(w'r'2)/dz
+IF (GFR2) THEN
+ ZF = ZF + M3_WTH_WR2(KKA,KKU,KKL,PREDTH1,PREDR1,PD,ZKEFF,PTKEM,&
+ & PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PEMOIST,PDTH_DZ) * MZM(PFR2, KKA, KKU, KKL)
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_WR2_O_DDTDZ(KKA,KKU,KKL,PREDTH1,PREDR1,PD,&
+ & ZKEFF,PTKEM,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PEMOIST) * MZM(PFR2, KKA, KKU, KKL)
+END IF
+!
+! d(w'th'r')/dz
+IF (GFTHR) THEN
+ Z3RDMOMENT= M3_WTH_WTHR(KKA,KKU,KKL,PREDR1,PD,ZKEFF,PTKEM,PSQRT_TKE,PBETA,&
+ & PLEPS,PEMOIST)
+!
+ ZF = ZF + Z3RDMOMENT * MZM(PFTHR, KKA, KKU, KKL)
+ ZDFDDTDZ = ZDFDDTDZ + D_M3_WTH_WTHR_O_DDTDZ(Z3RDMOMENT,PREDTH1,&
+ & PREDR1,PD,PBLL_O_E,PETHETA) * MZM(PFTHR, KKA, KKU, KKL)
+END IF
+!
+!* in 3DIM case, a part of the flux goes vertically, and another goes horizontally
+! (in presence of slopes)
+!* in 1DIM case, the part of energy released in horizontal flux
+! is taken into account in the vertical part
+!
+IF (HTURBDIM=='3DIM') THEN
+ ZF(:,:,IKB) = ( PIMPL*PSFTHP(:,:) + PEXPL*PSFTHM(:,:) ) &
+ * PDIRCOSZW(:,:) &
+ * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
+ELSE
+ ZF(:,:,IKB) = ( PIMPL*PSFTHP(:,:) + PEXPL*PSFTHM(:,:) ) &
+ / PDIRCOSZW(:,:) &
+ * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
+END IF
+!
+! Compute the splitted conservative potential temperature at t+deltat
+CALL TRIDIAG_THERMO(KKA,KKU,KKL,PTHLM,ZF,ZDFDDTDZ,PTSTEP,PIMPL,PDZZ,&
+ PRHODJ,PTHLP)
+!
+! Compute the equivalent tendency for the conservative potential temperature
+PRTHLS(:,:,:)= PRTHLS(:,:,:) + &
+ PRHODJ(:,:,:)*(PTHLP(:,:,:)-PTHLM(:,:,:))/PTSTEP
+!
+!* 2.2 Partial Thermal Production
+!
+! Conservative potential temperature flux :
+!
+ZFLXZ(:,:,:) = ZF &
+ + PIMPL * ZDFDDTDZ * DZM(PTHLP - PTHLM, KKA, KKU, KKL) / PDZZ
+!
+ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+!
+ DO JK=IKTB+1,IKTE-1
+ PWTH(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+KKL))
+ END DO
+ PWTH(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+KKL))
+ PWTH(:,:,KKA)=0.5*(ZFLXZ(:,:,KKA)+ZFLXZ(:,:,KKA+KKL))
+ PWTH(:,:,IKE)=PWTH(:,:,IKE-KKL)
+
+IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+ ! stores the conservative potential temperature vertical flux
+ TZFIELD%CMNHNAME = 'THW_FLX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'THW_FLX'
+ TZFIELD%CUNITS = 'K m s-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'Conservative potential temperature vertical flux'
+ TZFIELD%NGRID = 4
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
+END IF
+!
+! Contribution of the conservative temperature flux to the buoyancy flux
+IF (KRR /= 0) THEN
+ PTP(:,:,:) = PBETA * MZF(MZM(PETHETA, KKA, KKU, KKL) * ZFLXZ, KKA, KKU, KKL)
+ PTP(:,:,IKB)= PBETA(:,:,IKB) * PETHETA(:,:,IKB) * &
+ 0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
+ELSE
+ PTP(:,:,:)= PBETA * MZF(ZFLXZ, KKA, KKU, KKL)
+END IF
+!
+! Buoyancy flux at flux points
+!
+PWTHV = MZM(PETHETA, KKA, KKU, KKL) * ZFLXZ
+PWTHV(:,:,IKB) = PETHETA(:,:,IKB) * ZFLXZ(:,:,IKB)
+!
+!* 2.3 Partial vertical divergence of the < Rc w > flux
+! Correction for qc and qi negative in AROME
+!IF ( KRRL >= 1 ) THEN
+! IF ( KRRI >= 1 ) THEN
+! PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
+! DZF(MZM(PRHODJ*PATHETA*2.*PSRCM, KKA, KKU, KKL)*ZFLXZ/PDZZ, KKA, KKU, KKL) &
+! *(1.0-PFRAC_ICE(:,:,:))
+! PRRS(:,:,:,4) = PRRS(:,:,:,4) - &
+! DZF(MZM(PRHODJ*PATHETA*2.*PSRCM, KKA, KKU, KKL)*ZFLXZ/PDZZ, KKA, KKU, KKL) &
+! *PFRAC_ICE(:,:,:)
+! ELSE
+! PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
+! DZF(MZM(PRHODJ*PATHETA*2.*PSRCM, KKA, KKU, KKL)*ZFLXZ/PDZZ, KKA, KKU, KKL)
+! END IF
+!END IF
+!
+!* 2.4 Storage in LES configuration
+!
+IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WThl )
+ CALL LES_MEAN_SUBGRID(MZF(PWM*ZFLXZ, KKA, KKU, KKL), X_LES_RES_W_SBG_WThl )
+ CALL LES_MEAN_SUBGRID(GZ_W_M(PWM,PDZZ, KKA, KKU, KKL)*MZF(ZFLXZ, KKA, KKU, KKL),&
+ & X_LES_RES_ddxa_W_SBG_UaThl )
+ CALL LES_MEAN_SUBGRID(MZF(PDTH_DZ*ZFLXZ, KKA, KKU, KKL), X_LES_RES_ddxa_Thl_SBG_UaThl )
+ CALL LES_MEAN_SUBGRID(-XCTP*PSQRT_TKE/PLM*MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_ThlPz )
+ CALL LES_MEAN_SUBGRID(MZF(MZM(PETHETA, KKA, KKU, KKL)*ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WThv )
+ IF (KRR>=1) THEN
+ CALL LES_MEAN_SUBGRID(MZF(PDR_DZ*ZFLXZ, KKA, KKU, KKL), X_LES_RES_ddxa_Rt_SBG_UaThl )
+ END IF
+ !* diagnostic of mixing coefficient for heat
+ ZA = DZM(PTHLP, KKA, KKU, KKL)
+ WHERE (ZA==0.) ZA=1.E-6
+ ZA = - ZFLXZ / ZA * PDZZ
+ ZA(:,:,IKB) = XCSHF*PPHI3(:,:,IKB)*ZKEFF(:,:,IKB)
+ ZA = MZF(ZA, KKA, KKU, KKL)
+ ZA = MIN(MAX(ZA,-1000.),1000.)
+ CALL LES_MEAN_SUBGRID( ZA, X_LES_SUBGRID_Kh )
+ !
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+!
+!* 2.5 New boundary layer depth for TOMs
+!
+IF (HTOM=='TM06') CALL TM06_H(IKB,IKTB,IKTE,PTSTEP,PZZ,ZFLXZ,PBL_DEPTH)
+!
+!----------------------------------------------------------------------------
+!
+!
+!* 3. SOURCES OF CONSERVATIVE AND CLOUD MIXING RATIO AND
+! COMPLETE THERMAL PRODUCTION
+! ------------------------------------------------------
+!
+!* 3.1 Splitted value for cons. mixing ratio at t+deltat
+!
+!
+IF (KRR /= 0) THEN
+ ! Compute the turbulent flux F and F' at time t-dt.
+ !
+ IF (LHARAT) THEN
+ ZF (:,:,:) = -ZKEFF*DZM(PRM(:,:,:,1), KKA, KKU, KKL)/PDZZ
+ ZDFDDRDZ(:,:,:) = -ZKEFF
+ ELSE
+ ZF (:,:,:) = -XCSHF*PPSI3*ZKEFF*DZM(PRM(:,:,:,1), KKA, KKU, KKL)/PDZZ
+ ZDFDDRDZ(:,:,:) = -XCSHF*ZKEFF*D_PSI3DRDZ_O_DDRDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV)
+ ENDIF
+ !
+ ! Effect of 3rd order terms in temperature flux (at flux point)
+ !
+ ! d(w'2r')/dz
+ IF (GFWR) THEN
+ Z3RDMOMENT= M3_WR_W2R(KKA,KKU,KKL,PREDR1,PREDTH1,PD,ZKEFF,PTKEM)
+ !
+ ZF = ZF + Z3RDMOMENT * PFWR
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_W2R_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,PD,&
+ & PBLL_O_E,PEMOIST,ZKEFF,PTKEM) * PFWR
+ END IF
+ !
+ ! d(w'r'2)/dz
+ IF (GFR2) THEN
+ Z3RDMOMENT= M3_WR_WR2(PREDR1,PREDTH1,PD,PBLL_O_E,PEMOIST)
+ !
+ ZF = ZF + Z3RDMOMENT * MZM(PFR2, KKA, KKU, KKL)
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_WR2_O_DDRDZ(Z3RDMOMENT,PREDR1,&
+ & PREDTH1,PD,PBLL_O_E,PEMOIST) * MZM(PFR2, KKA, KKU, KKL)
+ END IF
+ !
+ ! d(w'2th')/dz
+ IF (GFWTH) THEN
+ ZF = ZF + M3_WR_W2TH(KKA,KKU,KKL,PREDR1,PREDTH1,PD,ZKEFF,&
+ & PTKEM,PBLL_O_E,PETHETA,PDR_DZ) * PFWTH
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_W2TH_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,&
+ & PD,ZKEFF,PTKEM,PBLL_O_E,PETHETA) * PFWTH
+ END IF
+ !
+ ! d(w'th'2)/dz
+ IF (GFTH2) THEN
+ ZF = ZF + M3_WR_WTH2(KKA,KKU,KKL,PREDR1,PREDTH1,PD,ZKEFF,PTKEM,&
+ & PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PETHETA,PDR_DZ) * MZM(PFTH2, KKA, KKU, KKL)
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_WTH2_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,PD,&
+ &ZKEFF,PTKEM,PSQRT_TKE,PBLL_O_E,PBETA,PLEPS,PETHETA) * MZM(PFTH2, KKA, KKU, KKL)
+ END IF
+ !
+ ! d(w'th'r')/dz
+ IF (GFTHR) THEN
+ Z3RDMOMENT= M3_WR_WTHR(KKA,KKU,KKL,PREDTH1,PD,ZKEFF,PTKEM,PSQRT_TKE,PBETA,&
+ & PLEPS,PETHETA)
+ !
+ ZF = ZF + Z3RDMOMENT * MZM(PFTHR, KKA, KKU, KKL)
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_WTHR_O_DDRDZ(KKA,KKU,KKL,Z3RDMOMENT,PREDR1, &
+ & PREDTH1,PD,PBLL_O_E,PEMOIST) * MZM(PFTHR, KKA, KKU, KKL)
+ END IF
+ !
+ !* in 3DIM case, a part of the flux goes vertically, and another goes horizontally
+ ! (in presence of slopes)
+ !* in 1DIM case, the part of energy released in horizontal flux
+ ! is taken into account in the vertical part
+ !
+ IF (HTURBDIM=='3DIM') THEN
+ ZF(:,:,IKB) = ( PIMPL*PSFRP(:,:) + PEXPL*PSFRM(:,:) ) &
+ * PDIRCOSZW(:,:) &
+ * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
+ ELSE
+ ZF(:,:,IKB) = ( PIMPL*PSFRP(:,:) + PEXPL*PSFRM(:,:) ) &
+ / PDIRCOSZW(:,:) &
+ * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
+ END IF
+ !
+ ! Compute the splitted conservative potential temperature at t+deltat
+ CALL TRIDIAG_THERMO(KKA,KKU,KKL,PRM(:,:,:,1),ZF,ZDFDDRDZ,PTSTEP,PIMPL,&
+ PDZZ,PRHODJ,PRP)
+ !
+ ! Compute the equivalent tendency for the conservative mixing ratio
+ PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRHODJ(:,:,:) * &
+ (PRP(:,:,:)-PRM(:,:,:,1))/PTSTEP
+ !
+ !* 3.2 Complete thermal production
+ !
+ ! cons. mixing ratio flux :
+ !
+ ZFLXZ(:,:,:) = ZF &
+ + PIMPL * ZDFDDRDZ * DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ
+ !
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ !
+ DO JK=IKTB+1,IKTE-1
+ PWRC(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+KKL))
+ END DO
+ PWRC(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+KKL))
+ PWRC(:,:,KKA)=0.5*(ZFLXZ(:,:,KKA)+ZFLXZ(:,:,KKA+KKL))
+ PWRC(:,:,IKE)=PWRC(:,:,IKE-KKL)
+ !
+ !
+ IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+ ! stores the conservative mixing ratio vertical flux
+ TZFIELD%CMNHNAME = 'RCONSW_FLX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'RCONSW_FLX'
+ TZFIELD%CUNITS = 'kg m s-1 kg-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'Conservative mixing ratio vertical flux'
+ TZFIELD%NGRID = 4
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
+ END IF
+ !
+ ! Contribution of the conservative water flux to the Buoyancy flux
+ ZA(:,:,:) = PBETA * MZF(MZM(PEMOIST, KKA, KKU, KKL) * ZFLXZ, KKA, KKU, KKL)
+ ZA(:,:,IKB) = PBETA(:,:,IKB) * PEMOIST(:,:,IKB) * &
+ 0.5 * ( ZFLXZ (:,:,IKB) + ZFLXZ (:,:,IKB+KKL) )
+ PTP(:,:,:) = PTP(:,:,:) + ZA(:,:,:)
+ !
+ ! Buoyancy flux at flux points
+ !
+ PWTHV = PWTHV + MZM(PEMOIST, KKA, KKU, KKL) * ZFLXZ
+ PWTHV(:,:,IKB) = PWTHV(:,:,IKB) + PEMOIST(:,:,IKB) * ZFLXZ(:,:,IKB)
+!
+!* 3.3 Complete vertical divergence of the < Rc w > flux
+! Correction of qc and qi negative for AROME
+! IF ( KRRL >= 1 ) THEN
+! IF ( KRRI >= 1 ) THEN
+! PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
+! DZF(MZM(PRHODJ*PAMOIST*2.*PSRCM, KKA, KKU, KKL)*ZFLXZ/PDZZ, KKA, KKU, KKL) &
+! *(1.0-PFRAC_ICE(:,:,:))
+! PRRS(:,:,:,4) = PRRS(:,:,:,4) - &
+! DZF(MZM(PRHODJ*PAMOIST*2.*PSRCM, KKA, KKU, KKL)*ZFLXZ/PDZZ, KKA, KKU, KKL) &
+! *PFRAC_ICE(:,:,:)
+! ELSE
+! PRRS(:,:,:,2) = PRRS(:,:,:,2) - &
+! DZF(MZM(PRHODJ*PAMOIST*2.*PSRCM, KKA, KKU, KKL)*ZFLXZ/PDZZ, KKA, KKU, KKL)
+! END IF
+! END IF
+!
+!* 3.4 Storage in LES configuration
+!
+ IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WRt )
+ CALL LES_MEAN_SUBGRID(MZF(PWM*ZFLXZ, KKA, KKU, KKL), X_LES_RES_W_SBG_WRt )
+ CALL LES_MEAN_SUBGRID(GZ_W_M(PWM,PDZZ, KKA, KKU, KKL)*MZF(ZFLXZ, KKA, KKU, KKL),&
+ & X_LES_RES_ddxa_W_SBG_UaRt )
+ CALL LES_MEAN_SUBGRID(MZF(PDTH_DZ*ZFLXZ, KKA, KKU, KKL), X_LES_RES_ddxa_Thl_SBG_UaRt )
+ CALL LES_MEAN_SUBGRID(MZF(PDR_DZ*ZFLXZ, KKA, KKU, KKL), X_LES_RES_ddxa_Rt_SBG_UaRt )
+ CALL LES_MEAN_SUBGRID(MZF(MZM(PEMOIST, KKA, KKU, KKL)*ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WThv , .TRUE. )
+ CALL LES_MEAN_SUBGRID(-XCTP*PSQRT_TKE/PLM*MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_RtPz )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+!
+END IF
+!
+!----------------------------------------------------------------------------
+!
+!
+!* 4. TURBULENT CORRELATIONS : <w Rc>
+! -------------------------------
+!
+!
+!* 4.1 <w Rc>
+!
+IF ( ((OTURB_FLX .AND. TPFILE%LOPENED) .OR. LLES_CALL) .AND. (KRRL > 0) ) THEN
+!
+! recover the Conservative potential temperature flux :
+! With LHARAT is true tke and length scales at half levels
+! yet modify to use length scale and tke at half levels from vdfexcuhl
+ IF (LHARAT) THEN
+ ZA(:,:,:) = DZM(PIMPL * PTHLP + PEXPL * PTHLM, KKA, KKU, KKL) / PDZZ * &
+ (-PLM*PSQRT_TKE)
+ ELSE
+ ZA(:,:,:) = DZM(PIMPL * PTHLP + PEXPL * PTHLM, KKA, KKU, KKL) / PDZZ * &
+ (-PPHI3*MZM(PLM*PSQRT_TKE, KKA, KKU, KKL)) * XCSHF
+ ENDIF
+ ZA(:,:,IKB) = ( PIMPL*PSFTHP(:,:) + PEXPL*PSFTHM(:,:) ) &
+ * PDIRCOSZW(:,:)
+ !
+ ! compute <w Rc>
+ ZFLXZ(:,:,:) = MZM(PAMOIST * 2.* PSRCM, KKA, KKU, KKL) * ZFLXZ(:,:,:) + &
+ MZM(PATHETA * 2.* PSRCM, KKA, KKU, KKL) * ZA(:,:,:)
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ !
+ ! store the liquid water mixing ratio vertical flux
+ IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+ TZFIELD%CMNHNAME = 'RCW_FLX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'RCW_FLX'
+ TZFIELD%CUNITS = 'kg m s-1 kg-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'Liquid water mixing ratio vertical flux'
+ TZFIELD%NGRID = 4
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
+ END IF
+ !
+! and we store in LES configuration this subgrid flux <w'rc'>
+!
+ IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WRc )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+!
+END IF !end of <w Rc>
+!
+!----------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_FLUX',1,ZHOOK_HANDLE)
+END SUBROUTINE TURB_VER_THERMO_FLUX
+END MODULE MODE_TURB_VER_THERMO_FLUX
--
GitLab