diff --git a/src/arome/aux/mode_argslist_ll.F90 b/src/arome/aux/mode_argslist_ll.F90
index ed5c4f1796c8fdb592c31a474a2b92c6b46d1b71..a01b2324990bdc736fc1884ed585fbf074677ecd 100644
--- a/src/arome/aux/mode_argslist_ll.F90
+++ b/src/arome/aux/mode_argslist_ll.F90
@@ -19,4 +19,15 @@ IMPLICIT NONE
!
CALL ABORT
END SUBROUTINE ADD3DFIELD_ll
+!
+ SUBROUTINE ADD4DFIELD_ll(TPLIST, PFIELD, HNAME)
+IMPLICIT NONE
+
+ TYPE(LIST_ll), POINTER :: TPLIST ! list of fields
+ REAL, DIMENSION(:,:,:,:), TARGET :: PFIELD ! field to be added to the list
+ ! of fields
+ character(len=*), intent(in) :: HNAME ! Name of the field to be added
+ !
+ CALL ABORT
+END SUBROUTINE ADD4DFIELD_ll
END MODULE MODE_ARGSLIST_ll
diff --git a/src/arome/aux/mode_ll.F90 b/src/arome/aux/mode_ll.F90
index a71daecfab08c9a546172b362a32210dec4b498c..58291ba77422036ecd9b200950125a028f023479 100644
--- a/src/arome/aux/mode_ll.F90
+++ b/src/arome/aux/mode_ll.F90
@@ -6,13 +6,12 @@ CONTAINS
SUBROUTINE GET_INDICE_ll(KXOR, KYOR, KXEND, KYEND, KSIZE1, KSIZE2)
USE MODD_PARAMETERS, ONLY : JPHEXT
IMPLICIT NONE
- INTEGER, INTENT(IN) :: KSIZE1, KSIZE2
+ INTEGER, INTENT(IN),OPTIONAL :: KSIZE1, KSIZE2
INTEGER, INTENT(OUT) :: KXOR, KYOR, KXEND, KYEND
KXOR=1+JPHEXT
KYOR=1+JPHEXT
KXEND=KSIZE1-JPHEXT
KYEND=KSIZE2-JPHEXT
- CALL ABORT
END SUBROUTINE GET_INDICE_ll
SUBROUTINE UPDATE_HALO_ll(TPLIST, KINFO)
diff --git a/src/arome/aux/mode_mppdb.F90 b/src/arome/aux/mode_mppdb.F90
new file mode 100644
index 0000000000000000000000000000000000000000..982b25d5deed69b423aa5ff6dae001cfe68099af
--- /dev/null
+++ b/src/arome/aux/mode_mppdb.F90
@@ -0,0 +1,18 @@
+MODULE MODE_MPPDB
+IMPLICIT NONE
+REAL :: PRECISION = 1e-8 * 0.0
+CONTAINS
+SUBROUTINE MPPDB_CHECK3DM(MESSAGE,PRECISION &
+ ,PTAB1,PTAB2,PTAB3,PTAB4,PTAB5,PTAB6,PTAB7,PTAB8,PTAB9,PTAB10 &
+ ,PTAB11,PTAB12,PTAB13,PTAB14,PTAB15,PTAB16,PTAB17,PTAB18,PTAB19,PTAB20 &
+ )
+
+IMPLICIT NONE
+
+CHARACTER(lEN=*) :: MESSAGE
+REAL :: PRECISION
+REAL, DIMENSION(:,:,:), OPTIONAL :: PTAB1,PTAB2,PTAB3,PTAB4,PTAB5,PTAB6,PTAB7,PTAB8,PTAB9,PTAB10
+REAL, DIMENSION(:,:,:), OPTIONAL :: PTAB11,PTAB12,PTAB13,PTAB14,PTAB15,PTAB16,PTAB17,PTAB18,PTAB19,PTAB20
+! DO NOTHING IN AROME
+END SUBROUTINE MPPDB_CHECK3DM
+END MODULE MODE_MPPDB
diff --git a/src/arome/aux/mode_tridiag_w.F90 b/src/arome/aux/mode_tridiag_w.F90
new file mode 100644
index 0000000000000000000000000000000000000000..6e5ce130493b6568b9888cbf865571da2c4421f2
--- /dev/null
+++ b/src/arome/aux/mode_tridiag_w.F90
@@ -0,0 +1,3 @@
+MODULE MODE_TRIDIAG_W
+! Empty module for PHYEX, used in TURB 3D
+END MODULE MODE_TRIDIAG_W
diff --git a/src/arome/aux/modi_shuman.F90 b/src/arome/aux/modi_shuman.F90
index 8bc69a410cc83ea136f24c444eb8abf328091418..d8ffd80a10bbe333b86fc74b43637d0024139fd7 100644
--- a/src/arome/aux/modi_shuman.F90
+++ b/src/arome/aux/modi_shuman.F90
@@ -35,16 +35,16 @@ END FUNCTION DYM
FUNCTION DZF(PA,KKA,KKU,KL) RESULT(PDZF)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux
! side
-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
+INTEGER, INTENT(IN),OPTIONAL :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN),OPTIONAL :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDZF ! result at mass localization
END FUNCTION DZF
!
FUNCTION DZM(PA,KKA,KKU,KL) RESULT(PDZM)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass
! localization
-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
+INTEGER, INTENT(IN),OPTIONAL :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN),OPTIONAL :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDZM ! result at flux side
END FUNCTION DZM
!
@@ -74,16 +74,16 @@ END FUNCTION MYM
!
FUNCTION MZF(PA,KKA,KKU,KL) RESULT(PMZF)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux side
-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
+INTEGER, INTENT(IN),OPTIONAL :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN),OPTIONAL :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZF ! result at mass
! localization
END FUNCTION MZF
!
FUNCTION MZM(PA,KKA,KKU,KL) RESULT(PMZM)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass localization
-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
+INTEGER, INTENT(IN),OPTIONAL :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN),OPTIONAL :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZM ! result at flux localization
END FUNCTION MZM
!
diff --git a/src/arome/aux/shuman.F90 b/src/arome/aux/shuman.F90
index 4aaa979958fc557c5501d0d06f1748012c23b79a..f8949e00d8e7966ffc2122154fc29121ec8f9a0e 100644
--- a/src/arome/aux/shuman.F90
+++ b/src/arome/aux/shuman.F90
@@ -88,7 +88,6 @@ PMXF=PA
!END DO
!
!PMXF(IIU,:,:) = PMXF(2*JPHEXT,:,:)
-CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
!
!-------------------------------------------------------------------------------
!
@@ -183,8 +182,6 @@ PMXM=PA
!
!PMXM(1,:,:) = PMXM(IIU-2*JPHEXT+1,:,:)
!
-CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
-
!-------------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('MXM',1,ZHOOK_HANDLE)
@@ -278,7 +275,6 @@ PMYF=PA
! PMYF(:,JJ,:) = 0.5*( PA(:,JJ,:)+PA(:,JJ+1,:) )
!END DO
!
-!PMYF(:,IJU,:) = PMYF(:,2*JPHEXT,:)
!
!-------------------------------------------------------------------------------
!
@@ -373,7 +369,6 @@ PMYM=PA
!
!PMYM(:,1,:) = PMYM(:,IJU-2*JPHEXT+1,:)
!
-CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
!-------------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('MYM',1,ZHOOK_HANDLE)
@@ -433,8 +428,8 @@ IMPLICIT NONE
! ------------------------------------
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux side
-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
+INTEGER, INTENT(IN),OPTIONAL :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN),OPTIONAL :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZF ! result at mass
! localization
!
@@ -517,8 +512,8 @@ IMPLICIT NONE
! ------------------------------------
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass localization
-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
+INTEGER, INTENT(IN),OPTIONAL :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN),OPTIONAL :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZM ! result at flux localization
!
!* 0.2 Declarations of local variables
@@ -626,7 +621,6 @@ DO JI=1,IIU-1
END DO
!
!PDXF(IIU,:,:) = PDXF(2*JPHEXT,:,:)
-CALL ABORT ! AROME SHOULD NOT CALLED HORIZONTAL FINITE DIFFERENCE
!
!-------------------------------------------------------------------------------
!
@@ -961,8 +955,8 @@ IMPLICIT NONE
! ------------------------------------
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at flux side
-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
+INTEGER, INTENT(IN),OPTIONAL :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN),OPTIONAL :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDZF ! result at mass
! localization
!
@@ -1045,8 +1039,8 @@ IMPLICIT NONE
! ------------------------------------
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! variable at mass localization
-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
+INTEGER, INTENT(IN),OPTIONAL :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER, INTENT(IN),OPTIONAL :: KL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDZM ! result at flux
! side
!
diff --git a/src/common/turb/mode_coefj.f90 b/src/common/turb/mode_coefj.f90
new file mode 100644
index 0000000000000000000000000000000000000000..239b3f63c3c5a3c774fd10014840d4e42bc9a17e
--- /dev/null
+++ b/src/common/turb/mode_coefj.f90
@@ -0,0 +1,135 @@
+!MNH_LIC Copyright 1994-2014 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.
+!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source$ $Revision$
+! MASDEV4_7 turb 2006/05/18 13:07:25
+!-----------------------------------------------------------------
+!################
+MODULE MODE_COEFJ
+IMPLICIT NONE
+CONTAINS
+! #######################################################
+ FUNCTION COEFJ(PTHL,PEXNREF,PFRAC_ICE) RESULT(PCOEFJ)
+! #######################################################
+!
+! PURPOSE
+!! -------
+! COEFJ computes the coefficient J of the documentation.
+!
+!!** METHOD
+!! ------ rvs(Tl) Lv(Tl)
+!! The value of this coefficient is J = --------------, for rc only
+!! Rv Tl THETAl
+!!
+!! rvsw(Tl) Lv(Tl) rvsi(Tl) Ls(Tl)
+!! or --------------- (1-Pfrac_ri) + --------------- Pfrac_ri, for rc+ri.
+!! Rv Tl THETAl Rv Tl THETAl
+!!
+!! EXTERNAL
+!! --------
+!! None.
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST : contains physical constants.
+!!
+!! REFERENCE
+!! ---------
+!! Book 1 of documentation of Meso-NH
+!! Book 2 of documentation of Meso-NH
+!!
+!!
+!! AUTHOR
+!! ------
+!! Jean-Marie Carriere * Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 20/03/95
+!! J.-P. Pinty 20/02/03 add non-precipitating ice
+!!
+!! ----------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of arguments and result
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHL ! Temperature variable
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Exner function of the
+! reference state
+REAL, DIMENSION(:,:,:), INTENT(IN), OPTIONAL :: PFRAC_ICE
+ ! Fraction of ri in the
+ ! non-precipating
+ ! "rc+ri" condensate
+!
+REAL,DIMENSION(SIZE(PTHL,1),SIZE(PTHL,2),SIZE(PTHL,3)):: PCOEFJ ! result
+!
+!* 0.2 declarations of local variables
+!
+REAL,DIMENSION(SIZE(PTHL,1),SIZE(PTHL,2),SIZE(PTHL,3)) :: &
+ ZTL, ZL, ZES, ZRVS, ZP
+! ZTL = Tl, ZL = Lv(Tl) or Ls(Tl), ZES = esw(Tl) or esi(Tl)
+! ZRVS = rvsw(Tl) or rvsi(Tl), ZP = p
+!
+REAL :: ZEPS ! = Mv/Md
+!---------------------------------------------------------------------------
+!
+!* 1. COMPUTATION OF Tl
+! -----------------
+!
+ZTL(:,:,:) = PTHL(:,:,:) * PEXNREF(:,:,:)
+!
+!* 2. COMPUTATION OF Lv(Tl)
+! ---------------------
+!
+ZL(:,:,:) = XLVTT + ( XCPV - XCL ) * ( ZTL(:,:,:) -XTT )
+!
+!* 3. COMPUTATION OF rvs(Tl)
+! ----------------------
+!
+ZEPS = XMV/XMD
+ZP(:,:,:) = (PEXNREF(:,:,:)**(XCPD/XRD))*XP00
+ZES(:,:,:) = EXP( XALPW - XBETAW/ZTL(:,:,:) - XGAMW*ALOG(ZTL(:,:,:) ) )
+ZRVS(:,:,:) = ZES(:,:,:) * ZEPS / ( ZP(:,:,:) - ZES(:,:,:) )
+!
+! 4. RESULT FOR rc only
+! ------------------
+!
+PCOEFJ(:,:,:) = ZRVS(:,:,:)*ZL(:,:,:)/ &
+ ( XRV*ZTL(:,:,:)*PTHL(:,:,:) )
+!
+! Add case when rc+ri
+!
+IF(PRESENT(PFRAC_ICE)) THEN
+!
+!* 5. COMPUTATION OF Ls(Tl)
+! ---------------------
+!
+ ZL(:,:,:) = XLSTT + ( XCPV - XCI ) * ( ZTL(:,:,:) -XTT )
+!
+!* 6. COMPUTATION OF rvs(Tl)
+! ----------------------
+!
+ ZES(:,:,:) = EXP( XALPI - XBETAI/ZTL(:,:,:) - XGAMI*ALOG(ZTL(:,:,:) ) )
+ ZRVS(:,:,:) = ZES(:,:,:) * ZEPS / ( ZP(:,:,:) - ZES(:,:,:) )
+!
+! 7. RESULT FOR rc and ri
+! --------------------
+!
+ PCOEFJ(:,:,:) = (1.0 - PFRAC_ICE(:,:,:))*PCOEFJ(:,:,:) &
+ + PFRAC_ICE(:,:,:) *ZRVS(:,:,:)*ZL(:,:,:)/ &
+ ( XRV*ZTL(:,:,:)*PTHL(:,:,:) )
+END IF
+!
+!---------------------------------------------------------------------------
+!
+END FUNCTION COEFJ
+END MODULE MODE_COEFJ
diff --git a/src/common/turb/mode_turb_hor.F90 b/src/common/turb/mode_turb_hor.F90
new file mode 100644
index 0000000000000000000000000000000000000000..993b60a5539e87e6b164f9d35a9b9e1f73858880
--- /dev/null
+++ b/src/common/turb/mode_turb_hor.F90
@@ -0,0 +1,367 @@
+!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_HOR
+IMPLICIT NONE
+CONTAINS
+ SUBROUTINE TURB_HOR(KSPLT, KRR, KRRL, KRRI, PTSTEP, &
+ OTURB_FLX,OSUBG_COND, &
+ TPFILE, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, &
+ PDIRCOSXW,PDIRCOSYW,PDIRCOSZW, &
+ PCOSSLOPE,PSINSLOPE, &
+ PINV_PDXX, PINV_PDYY, PINV_PDZZ, PMZM_PRHODJ, &
+ PK, &
+ PRHODJ,PTHVREF, &
+ PSFTHM,PSFRM,PSFSVM, &
+ PCDUEFF,PTAU11M,PTAU12M,PTAU22M,PTAU33M, &
+ PUM,PVM,PWM,PUSLOPEM,PVSLOPEM,PTHLM,PRM,PSVM, &
+ PTKEM,PLM,PLEPS, &
+ PLOCPEXNM,PATHETA,PAMOIST,PSRCM,PFRAC_ICE, &
+ PDP,PTP,PSIGS, &
+ PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS )
+! ################################################################
+!
+!
+!!**** *TURB_HOR* -routine to compute the source terms in the meso-NH
+!! model equations due to the non-vertical turbulent fluxes.
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute the non-vertical
+! turbulent fluxes of the evolutive variables and give back the
+! source terms to the main program.
+!
+!!** METHOD
+!! ------
+!! Complementary 3D calculations when running at high resolution;
+!! The non-vertical turbulent fluxes are computed explicitly. The
+!! contributions are cumulated in PRvarS and in DP and TP of TKE
+!
+! d(rho*T) = -d(rho*u'T'/dxx) -d(-rho*u'T'*dzx/dxx/dzz)
+! / dt / dx /dz
+!!
+!!
+!! Near the bottom of the model, uncentred evaluation of vertical
+!! gradients are required because no field values are available under
+!! the level where the gradient must be evaluated. In this case, the
+!! gradient is computed with a second order accurate uncentred scheme
+!! according to:
+!!
+!! D FF dzz3 (dzz3+dzz4)
+!! ---- = - ----------------- FF(4) + ----------------- FF(3)
+!! D z (dzz3+dzz4) dzz4 dzz3 dzz4
+!!
+!! dzz4 + 2 dzz3
+!! - ----------------- FF(2)
+!! (dzz3+dzz4) dzz3
+!!
+!! where the values are taken from:
+!!
+!! ----- FF(5)
+!! |
+!! | dzz5
+!! |
+!! ----- FF(4)
+!! |
+!! | dzz4
+!! |
+!! ----- FF(3)
+!! |
+!! | dzz3
+!! |
+!! ----- FF(2) , (D FF / DZ)
+!! | dzz2 * 0.5
+!! ----- ground
+!!
+!!
+!!
+!! EXTERNAL
+!! --------
+!!
+!!
+!! 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 : number of vertical external points
+!!
+!!
+!!
+!!
+!! REFERENCE
+!! ---------
+!! Book 2 of documentation (routine TURB_HOR)
+!! Book 1 of documentation (Chapter: Turbulence)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original Aug 29, 1994
+!! Modifications: Feb 14, 1995 (J.Cuxart and J.Stein)
+!! Doctorization and Optimization
+!! March 21, 1995 (J.M. Carriere)
+!! Introduction of cloud water
+!! June 14, 1995 (J. Stein)
+!! rm the ZVTPV computation + bug in the all
+!! or nothing condens. case
+!! June 28, 1995 (J.Cuxart) Add the LES tools
+!! Sept 19, 1995 (J. Stein) change the surface flux
+!! computations
+!! Nov 13, 1995 (J. Stein) include the tangential fluxes
+!! bug in <u'w'> at the surface
+!! Nov 27, 1997 (V. Saravane) spliting of the routine
+!! Nov 27, 1997 (V. Masson) clearing of the routine
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! Feb 20, 2003 (JP Pinty) Add PFRAC_ICE
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+USE MODD_CTURB
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_PARAMETERS
+USE MODD_LES
+!
+USE MODE_TURB_HOR_THERMO_FLUX, ONLY: TURB_HOR_THERMO_FLUX
+USE MODE_TURB_HOR_THERMO_CORR, ONLY: TURB_HOR_THERMO_CORR
+USE MODE_TURB_HOR_DYN_CORR, ONLY: TURB_HOR_DYN_CORR
+USE MODE_TURB_HOR_UV, ONLY: TURB_HOR_UV
+USE MODE_TURB_HOR_UW, ONLY: TURB_HOR_UW
+USE MODE_TURB_HOR_VW, ONLY: TURB_HOR_VW
+USE MODE_TURB_HOR_SV_FLUX, ONLY: TURB_HOR_SV_FLUX
+USE MODE_TURB_HOR_SV_CORR, ONLY: TURB_HOR_SV_CORR
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declaration of arguments
+!
+!
+INTEGER, INTENT(IN) :: KSPLT ! current split index
+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.
+REAL, INTENT(IN) :: PTSTEP !
+LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
+ ! turbulent fluxes in the syncronous FM-file
+LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for sub-grid
+! condensation
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! vertical grid
+REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSXW, PDIRCOSYW, PDIRCOSZW
+! Director Cosinus along x, y and z directions at surface w-point
+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) :: PK ! Turbulent diffusion doef.
+ ! PK = PLM * SQRT(PTKEM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDXX ! 1./PDXX
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDYY ! 1./PDYY
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDZZ ! 1./PDZZ
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PMZM_PRHODJ ! MZM(PRHODJ)
+
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density * grid volume
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! ref. state VPT
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHM,PSFRM
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSFSVM ! surface fluxes
+!
+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) :: PTAU22M ! <vv> in the same axes
+REAL, DIMENSION(:,:), INTENT(IN) :: PTAU33M ! <ww> in the same axes
+!
+! Variables at t-1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUM,PVM,PWM,PTHLM
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! mixing ratios at t-1,
+ ! where PRM(:,:,:,1) = conservative mixing ratio
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! scalar var. at t-1
+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- dt
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLEPS ! dissipative length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLOCPEXNM ! Lv(T)/Cp/Exner 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(INOUT) :: PRUS, PRVS, PRWS, PRTHLS
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRSVS,PRRS ! var. at t+1 -split-
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDP,PTP ! TKE production terms
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSIGS
+ ! IN: Vertical part of Sigma_s at t
+ ! OUT: Total Sigma_s at t
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+! ---------------------------------------------------------------------------
+!
+!* 1. PRELIMINARY COMPUTATIONS
+! ------------------------
+!
+!* Exchange coefficient is limited in order to insure numerical stability
+!
+!!
+!* 2. < U' THETA'l >
+!* 3. < U' R'np >
+!* 4. < U' TPV' >
+!* 5. < V' THETA'l >
+!* 6. < V' R'np >
+!* 7. < V' TPV' >
+!
+ CALL TURB_HOR_THERMO_FLUX(KSPLT, KRR, KRRL, KRRI, &
+ OTURB_FLX,OSUBG_COND, &
+ TPFILE, &
+ PK,PINV_PDXX,PINV_PDYY,PINV_PDZZ,PMZM_PRHODJ, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ PDIRCOSXW,PDIRCOSYW, &
+ PRHODJ, &
+ PSFTHM,PSFRM, &
+ PWM,PTHLM,PRM, &
+ PATHETA,PAMOIST,PSRCM,PFRAC_ICE, &
+ PRTHLS,PRRS )
+!
+!
+!* 8. TURBULENT CORRELATIONS : <THl THl>, <THl Rnp>, <Rnp Rnp>, Sigma_s
+!
+ IF (KSPLT==1) &
+ CALL TURB_HOR_THERMO_CORR(KRR, KRRL, KRRI, &
+ OTURB_FLX,OSUBG_COND, &
+ TPFILE, &
+ PINV_PDXX,PINV_PDYY, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ PTHVREF, &
+ PWM,PTHLM,PRM, &
+ PTKEM,PLM,PLEPS, &
+ PLOCPEXNM,PATHETA,PAMOIST,PSRCM, &
+ PSIGS )
+!
+!
+!* 9. < U'U'>
+!* 10. < V'V'>
+!* 11. < W'W'>
+!
+ CALL TURB_HOR_DYN_CORR(KSPLT, PTSTEP, &
+ OTURB_FLX,KRR, &
+ TPFILE, &
+ PK,PINV_PDZZ, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, &
+ PDIRCOSZW, &
+ PCOSSLOPE,PSINSLOPE, &
+ PRHODJ, &
+ PCDUEFF,PTAU11M,PTAU12M,PTAU22M,PTAU33M, &
+ PUM,PVM,PWM, PUSLOPEM,PVSLOPEM, &
+ PTHLM,PRM,PSVM, &
+ PTKEM,PLM, &
+ PDP,PTP, &
+ PRUS,PRVS,PRWS )
+!
+!
+!* 12. < U'V'>
+!
+ CALL TURB_HOR_UV(KSPLT, &
+ OTURB_FLX, &
+ TPFILE, &
+ PK,PINV_PDXX,PINV_PDYY,PINV_PDZZ,PMZM_PRHODJ, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ PDIRCOSZW, &
+ PCOSSLOPE,PSINSLOPE, &
+ PRHODJ, &
+ PCDUEFF,PTAU11M,PTAU12M,PTAU22M,PTAU33M, &
+ PUM,PVM,PUSLOPEM,PVSLOPEM, &
+ PDP, &
+ PRUS,PRVS )
+!
+!
+!* 13. < U'W'>
+!
+ CALL TURB_HOR_UW(KSPLT, &
+ OTURB_FLX,KRR, &
+ TPFILE, &
+ PK,PINV_PDXX,PINV_PDZZ,PMZM_PRHODJ, &
+ PDXX,PDZZ,PDZX, &
+ PRHODJ,PTHVREF, &
+ PUM,PWM,PTHLM,PRM,PSVM, &
+ PTKEM,PLM, &
+ PDP, &
+ PRUS,PRWS )
+!
+!
+!* 14. < V'W'>
+!
+ CALL TURB_HOR_VW(KSPLT, &
+ OTURB_FLX,KRR, &
+ TPFILE, &
+ PK,PINV_PDYY,PINV_PDZZ,PMZM_PRHODJ, &
+ PDYY,PDZZ,PDZY, &
+ PRHODJ,PTHVREF, &
+ PVM,PWM,PTHLM,PRM,PSVM, &
+ PTKEM,PLM, &
+ PDP, &
+ PRVS,PRWS )
+
+!
+!
+!* 15. HORIZONTAL FLUXES OF PASSIVE SCALARS
+!
+ CALL TURB_HOR_SV_FLUX(KSPLT, &
+ OTURB_FLX, &
+ TPFILE, &
+ PK,PINV_PDXX,PINV_PDYY,PINV_PDZZ,PMZM_PRHODJ, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ PDIRCOSXW,PDIRCOSYW, &
+ PRHODJ,PWM, &
+ PSFSVM, &
+ PSVM, &
+ PRSVS )
+!
+ IF (KSPLT==1 .AND. LLES_CALL) &
+ CALL TURB_HOR_SV_CORR(KRR,KRRL,KRRI, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ PLM,PLEPS,PTKEM,PTHVREF, &
+ PTHLM,PRM, &
+ PLOCPEXNM,PATHETA,PAMOIST,PSRCM, &
+ PWM,PSVM )
+!
+!
+END SUBROUTINE TURB_HOR
+END MODULE MODE_TURB_HOR
diff --git a/src/common/turb/mode_turb_hor_dyn_corr.F90 b/src/common/turb/mode_turb_hor_dyn_corr.F90
new file mode 100644
index 0000000000000000000000000000000000000000..674f91bdca45afb3f4dd12ef3b111009471491bf
--- /dev/null
+++ b/src/common/turb/mode_turb_hor_dyn_corr.F90
@@ -0,0 +1,556 @@
+!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_HOR_DYN_CORR
+IMPLICIT NONE
+CONTAINS
+ SUBROUTINE TURB_HOR_DYN_CORR(KSPLT, PTSTEP, &
+ OTURB_FLX,KRR, &
+ TPFILE, &
+ PK,PINV_PDZZ, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, &
+ PDIRCOSZW, &
+ PCOSSLOPE,PSINSLOPE, &
+ PRHODJ, &
+ PCDUEFF,PTAU11M,PTAU12M,PTAU22M,PTAU33M, &
+ PUM,PVM,PWM,PUSLOPEM,PVSLOPEM, &
+ PTHLM,PRM,PSVM, &
+ PTKEM,PLM, &
+ PDP,PTP, &
+ PRUS,PRVS,PRWS )
+! ################################################################
+!
+!!**** *TURB_HOR* -routine to compute the source terms in the meso-NH
+!! model equations due to the non-vertical turbulent fluxes.
+!!
+!! PURPOSE
+!! -------
+!!
+!! see TURB_HOR
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!!
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Aug , 1997 (V. Saravane) spliting of TURB_HOR
+!! Nov 27, 1997 (V. Masson) clearing of the routine
+!! Oct 18, 2000 (V. Masson) LES computations + LFLAT switch
+!! Feb 15, 2001 (J. Stein) remove the use of w=0 at the
+!! ground
+!! Mar 12, 2001 (V. Masson and J. Stein) major bugs
+!! + change of discretization at the surface
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
+!! change of YCOMMENT
+!! July 2012 (V.Masson) Implicitness of W
+!! March 2014 (V.Masson) tridiag_w : bug between
+!! mass and flux position
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_ARGSLIST_ll, ONLY: LIST_ll
+USE MODD_CST
+USE MODD_CONF
+USE MODD_CTURB
+USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_PARAMETERS
+USE MODD_LES
+USE MODD_NSV
+!
+USE MODE_ll
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
+!
+USE MODI_GRADIENT_M
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_SHUMAN
+USE MODE_COEFJ, ONLY: COEFJ
+USE MODI_LES_MEAN_SUBGRID
+USE MODE_TRIDIAG_W
+!
+USE MODI_SECOND_MNH
+USE MODE_MPPDB
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declaration of arguments
+!
+!
+!
+INTEGER, INTENT(IN) :: KSPLT ! split process index
+REAL, INTENT(IN) :: PTSTEP ! timestep
+LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
+ ! turbulent fluxes in the syncronous FM-file
+INTEGER, INTENT(IN) :: KRR ! number of moist var.
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PK ! Turbulent diffusion doef.
+ ! PK = PLM * SQRT(PTKEM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDZZ ! 1./PDZZ
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! vertical grid
+REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSZW
+! Director Cosinus along z directions at surface w-point
+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 ! density * grid volume
+!
+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) :: PTAU22M ! <vv> in the same axes
+REAL, DIMENSION(:,:), INTENT(IN) :: PTAU33M ! <ww> in the same axes
+!
+! Variables at t-1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUM,PVM,PWM,PTHLM
+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- dt
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS, PRWS
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDP,PTP ! TKE production terms
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),SIZE(PUM,3)) &
+ :: ZFLX,ZWORK
+ ! work arrays, PK is the turb. mixing coef.
+!
+REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2)) ::ZDIRSINZW
+ ! sinus of the angle between the vertical and the normal to the orography
+INTEGER :: IKB,IKE
+ ! Index values for the Beginning and End
+ ! mass points of the domain
+INTEGER :: IKU
+INTEGER :: JSV ! scalar loop counter
+!
+REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),SIZE(PUM,3)) :: GX_U_M_PUM
+REAL, DIMENSION(SIZE(PVM,1),SIZE(PVM,2),SIZE(PVM,3)) :: GY_V_M_PVM
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),SIZE(PWM,3)) :: GZ_W_M_PWM
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),SIZE(PWM,3)) :: GZ_W_M_ZWP
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),SIZE(PWM,3)) :: ZMZF_DZZ ! MZF(PDZZ)
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),SIZE(PWM,3)) :: ZDFDDWDZ ! formal derivative of the
+! ! flux (variable: dW/dz)
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),SIZE(PWM,3)) :: ZWP ! W at future time-step
+!
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),1) :: ZDU_DZ_DZS_DX ! du/dz*dzs/dx surf
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),1) :: ZDV_DZ_DZS_DY ! dv/dz*dzs/dy surf
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),1) :: ZDU_DX ! du/dx surf
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),1) :: ZDV_DY ! dv/dy surf
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),1) :: ZDW_DZ ! dw/dz surf
+!
+INTEGER :: IINFO_ll ! return code of parallel routine
+TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange
+
+REAL :: ZTIME1, ZTIME2
+
+
+REAL, DIMENSION(SIZE(PDZZ,1),SIZE(PDZZ,2),1+JPVEXT:3+JPVEXT) :: ZCOEFF , ZDZZ
+ ! coefficients for the uncentred gradient
+ ! computation near the ground
+TYPE(TFIELDDATA) :: TZFIELD
+! --------------------------------------------------------------------------
+!
+!* 1. PRELIMINARY COMPUTATIONS
+! ------------------------
+NULLIFY(TZFIELDS_ll)
+!
+IKB = 1+JPVEXT
+IKE = SIZE(PUM,3)-JPVEXT
+IKU = SIZE(PUM,3)
+!
+!
+ZDIRSINZW(:,:) = SQRT( 1. - PDIRCOSZW(:,:)**2 )
+!
+GX_U_M_PUM = GX_U_M(PUM,PDXX,PDZZ,PDZX)
+IF (.NOT. L2D) GY_V_M_PVM = GY_V_M(PVM,PDYY,PDZZ,PDZY)
+GZ_W_M_PWM = GZ_W_M(PWM,PDZZ)
+!
+ZMZF_DZZ = MZF(PDZZ)
+!
+CALL ADD3DFIELD_ll( TZFIELDS_ll, ZFLX, 'TURB_HOR_DYN_CORR::ZFLX' )
+
+
+! compute the coefficients for the uncentred gradient computation near the
+! ground
+!
+!* 9. < U'U'>
+! -------
+!
+! Computes the U variance
+IF (.NOT. L2D) THEN
+ ZFLX(:,:,:)= (2./3.) * PTKEM &
+ - XCMFS * PK *( (4./3.) * GX_U_M_PUM &
+ -(2./3.) * ( GY_V_M_PVM &
+ +GZ_W_M_PWM ) )
+ !! & to be tested later
+ !! + XCMFB * PLM / SQRT(PTKEM) * (-2./3.) * PTP
+ELSE
+ ZFLX(:,:,:)= (2./3.) * PTKEM &
+ - XCMFS * PK *( (4./3.) * GX_U_M_PUM &
+ -(2./3.) * ( GZ_W_M_PWM ) )
+ !! & to be tested later
+ !! + XCMFB * PLM / SQRT(PTKEM) * (-2./3.) * PTP
+END IF
+!
+ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
+!
+!* prescription of du/dz and dv/dz with uncentered gradient at the surface
+! prescription of dw/dz at Dz/2 above ground using the continuity equation
+! using a Boussinesq hypothesis to remove the z dependance of rhod_ref
+! (div u = 0)
+!
+ZDZZ(:,:,:) = MXM(PDZZ(:,:,IKB:IKB+2))
+ZCOEFF(:,:,IKB+2)= - ZDZZ(:,:,2) / &
+ ( (ZDZZ(:,:,3)+ZDZZ(:,:,2)) * ZDZZ(:,:,3) )
+ZCOEFF(:,:,IKB+1)= (ZDZZ(:,:,3)+ZDZZ(:,:,2)) / &
+ ( ZDZZ(:,:,2) * ZDZZ(:,:,3) )
+ZCOEFF(:,:,IKB)= - (ZDZZ(:,:,3)+2.*ZDZZ(:,:,2)) / &
+ ( (ZDZZ(:,:,3)+ZDZZ(:,:,2)) * ZDZZ(:,:,2) )
+!
+ZDU_DZ_DZS_DX(:,:,:)=MXF ((ZCOEFF(:,:,IKB+2:IKB+2)*PUM(:,:,IKB+2:IKB+2) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PUM(:,:,IKB+1:IKB+1) &
+ +ZCOEFF(:,:,IKB :IKB )*PUM(:,:,IKB :IKB ) &
+ )* 0.5 * ( PDZX(:,:,IKB+1:IKB+1)+PDZX(:,:,IKB:IKB)) &
+ )/ MXF(PDXX(:,:,IKB:IKB))
+!
+ZDZZ(:,:,:) = MYM(PDZZ(:,:,IKB:IKB+2))
+ZCOEFF(:,:,IKB+2)= - ZDZZ(:,:,2) / &
+ ( (ZDZZ(:,:,3)+ZDZZ(:,:,2)) * ZDZZ(:,:,3) )
+ZCOEFF(:,:,IKB+1)= (ZDZZ(:,:,3)+ZDZZ(:,:,2)) / &
+ ( ZDZZ(:,:,2) * ZDZZ(:,:,3) )
+ZCOEFF(:,:,IKB)= - (ZDZZ(:,:,3)+2.*ZDZZ(:,:,2)) / &
+ ( (ZDZZ(:,:,3)+ZDZZ(:,:,2)) * ZDZZ(:,:,2) )
+!
+
+ZDV_DZ_DZS_DY(:,:,:)=MYF ((ZCOEFF(:,:,IKB+2:IKB+2)*PVM(:,:,IKB+2:IKB+2) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PVM(:,:,IKB+1:IKB+1) &
+ +ZCOEFF(:,:,IKB :IKB )*PVM(:,:,IKB :IKB ) &
+ )* 0.5 * ( PDZY(:,:,IKB+1:IKB+1)+PDZY(:,:,IKB:IKB)) &
+ )/ MYF(PDYY(:,:,IKB:IKB))
+!
+!
+ZDU_DX(:,:,:)= DXF(PUM(:,:,IKB:IKB)) / MXF(PDXX(:,:,IKB:IKB)) &
+ - ZDU_DZ_DZS_DX(:,:,:)
+
+ZDV_DY(:,:,:)= DYF(PVM(:,:,IKB:IKB)) / MYF(PDYY(:,:,IKB:IKB)) &
+ - ZDV_DZ_DZS_DY(:,:,:)
+!
+ZDW_DZ(:,:,:)=-ZDU_DX(:,:,:)-ZDV_DY(:,:,:)
+!
+!* computation
+!
+ZFLX(:,:,IKB) = (2./3.) * PTKEM(:,:,IKB) &
+ - XCMFS * PK(:,:,IKB) * 2. * ZDU_DX(:,:,1)
+
+
+!! & to be tested later
+!! + XCMFB * PLM(:,:,IKB:IKB) /SQRT(PTKEM(:,:,IKB:IKB)) * &
+!! (-2./3.) * PTP(:,:,IKB:IKB)
+!
+! extrapolates this flux under the ground with the surface flux
+ZFLX(:,:,IKB-1) = &
+ PTAU11M(:,:) * PCOSSLOPE(:,:)**2 * PDIRCOSZW(:,:)**2 &
+ -2. * PTAU12M(:,:) * PCOSSLOPE(:,:)* PSINSLOPE(:,:) * PDIRCOSZW(:,:) &
+ + PTAU22M(:,:) * PSINSLOPE(:,:)**2 &
+ + PTAU33M(:,:) * PCOSSLOPE(:,:)**2 * ZDIRSINZW(:,:)**2 &
+ +2. * PCDUEFF(:,:) * ( &
+ PVSLOPEM(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) &
+ - PUSLOPEM(:,:) * PCOSSLOPE(:,:)**2 * ZDIRSINZW(:,:) * PDIRCOSZW(:,:) )
+!
+ZFLX(:,:,IKB-1) = 2. * ZFLX(:,:,IKB-1) - ZFLX(:,:,IKB)
+!
+CALL UPDATE_HALO_ll(TZFIELDS_ll, IINFO_ll)
+IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+ ! stores <U U>
+ TZFIELD%CMNHNAME = 'U_VAR'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'U_VAR'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_U_VAR'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+END IF
+!
+! Complete the U tendency
+IF (.NOT. LFLAT) THEN
+CALL MPPDB_CHECK3DM("before turb_corr:PRUS,PRHODJ,ZFLX,PDXX,PDZX,PINV_PDZZ",PRECISION,&
+ & PRUS,PRHODJ,ZFLX,PDXX,PDZX,PINV_PDZZ )
+
+ PRUS(:,:,:)=PRUS &
+ -DXM(PRHODJ * ZFLX / MXF(PDXX) ) &
+ +DZF( PDZX / MZM(PDXX) * MXM( MZM(PRHODJ*ZFLX) * PINV_PDZZ ) )
+CALL MPPDB_CHECK3DM("after turb_corr:PRUS,PRHODJ,ZFLX,PDXX,PDZX,PINV_PDZZ",PRECISION,&
+ & PRUS,PRHODJ,ZFLX,PDXX,PDZX,PINV_PDZZ )
+ELSE
+ PRUS(:,:,:)=PRUS -DXM(PRHODJ * ZFLX / MXF(PDXX) )
+END IF
+!
+IF (KSPLT==1) THEN
+ ! Contribution to the dynamic production of TKE:
+ ZWORK(:,:,:) = - ZFLX(:,:,:) * GX_U_M_PUM
+ !
+ ! evaluate the dynamic production at w(IKB+1) in PDP(IKB)
+ !
+ ZWORK(:,:,IKB) = 0.5* ( -ZFLX(:,:,IKB)*ZDU_DX(:,:,1) + ZWORK(:,:,IKB+1) )
+ !
+ PDP(:,:,:) = PDP(:,:,:) + ZWORK(:,:,:)
+END IF
+!
+! Storage in the LES configuration
+!
+IF (LLES_CALL .AND. KSPLT==1) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( ZFLX, X_LES_SUBGRID_U2 )
+ CALL LES_MEAN_SUBGRID( -ZWORK, X_LES_RES_ddxa_U_SBG_UaU , .TRUE.)
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+
+!
+!* 10. < V'V'>
+! -------
+!
+! Computes the V variance
+IF (.NOT. L2D) THEN
+ ZFLX(:,:,:)= (2./3.) * PTKEM &
+ - XCMFS * PK *( (4./3.) * GY_V_M_PVM &
+ -(2./3.) * ( GX_U_M_PUM &
+ +GZ_W_M_PWM ) )
+ !! & to be tested
+ !! + XCMFB * PLM / SQRT(PTKEM) * (-2./3.) * PTP
+ !
+ELSE
+ ZFLX(:,:,:)= (2./3.) * PTKEM &
+ - XCMFS * PK *(-(2./3.) * ( GX_U_M_PUM &
+ +GZ_W_M_PWM ) )
+ !! & to be tested
+ !! + XCMFB * PLM / SQRT(PTKEM) * (-2./3.) * PTP
+ !
+END IF
+!
+ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
+!
+ZFLX(:,:,IKB) = (2./3.) * PTKEM(:,:,IKB) &
+ - XCMFS * PK(:,:,IKB) * 2. * ZDV_DY(:,:,1)
+
+!! & to be tested
+!! + XCMFB * PLM(:,:,IKB:IKB) /SQRT(PTKEM(:,:,IKB:IKB)) * &
+!! (-2./3.) * PTP(:,:,IKB:IKB)
+!
+! extrapolates this flux under the ground with the surface flux
+ZFLX(:,:,IKB-1) = &
+ PTAU11M(:,:) * PSINSLOPE(:,:)**2 * PDIRCOSZW(:,:)**2 &
+ +2. * PTAU12M(:,:) * PCOSSLOPE(:,:)* PSINSLOPE(:,:) * PDIRCOSZW(:,:) &
+ + PTAU22M(:,:) * PCOSSLOPE(:,:)**2 &
+ + PTAU33M(:,:) * PSINSLOPE(:,:)**2 * ZDIRSINZW(:,:)**2 &
+ -2. * PCDUEFF(:,:)* ( &
+ PUSLOPEM(:,:) * PSINSLOPE(:,:)**2 * ZDIRSINZW(:,:) * PDIRCOSZW(:,:) &
+ + PVSLOPEM(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) )
+!
+ZFLX(:,:,IKB-1) = 2. * ZFLX(:,:,IKB-1) - ZFLX(:,:,IKB)
+!
+CALL UPDATE_HALO_ll(TZFIELDS_ll, IINFO_ll)
+!
+IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+ ! stores <V V>
+ TZFIELD%CMNHNAME = 'V_VAR'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'V_VAR'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_V_VAR'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+END IF
+!
+! Complete the V tendency
+IF (.NOT. L2D) THEN
+ IF (.NOT. LFLAT) THEN
+ PRVS(:,:,:)=PRVS &
+ -DYM(PRHODJ * ZFLX / MYF(PDYY) ) &
+ +DZF( PDZY / MZM(PDYY) * &
+ MYM( MZM(PRHODJ*ZFLX) * PINV_PDZZ ) )
+ ELSE
+ PRVS(:,:,:)=PRVS -DYM(PRHODJ * ZFLX / MYF(PDYY) )
+ END IF
+!
+! Contribution to the dynamic production of TKE:
+ IF (KSPLT==1) ZWORK(:,:,:) = - ZFLX(:,:,:) * GY_V_M_PVM
+ELSE
+ ZWORK(:,:,:) = 0.
+END IF
+!
+IF (KSPLT==1) THEN
+ !
+ ! evaluate the dynamic production at w(IKB+1) in PDP(IKB)
+ !
+ ZWORK(:,:,IKB) = 0.5* ( -ZFLX(:,:,IKB)*ZDV_DY(:,:,1) + ZWORK(:,:,IKB+1) )
+ !
+ PDP(:,:,:) = PDP(:,:,:) + ZWORK(:,:,:)
+END IF
+!
+! Storage in the LES configuration
+!
+IF (LLES_CALL .AND. KSPLT==1) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( ZFLX, X_LES_SUBGRID_V2 )
+ CALL LES_MEAN_SUBGRID( -ZWORK, X_LES_RES_ddxa_V_SBG_UaV , .TRUE.)
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+!
+!* 11. < W'W'>
+! -------
+!
+! Computes the W variance
+IF (.NOT. L2D) THEN
+ ZFLX(:,:,:)= (2./3.) * PTKEM &
+ - XCMFS * PK *( (4./3.) * GZ_W_M_PWM &
+ -(2./3.) * ( GX_U_M_PUM &
+ +GY_V_M_PVM ) )
+ !! & to be tested
+ !! -2.* XCMFB * PLM / SQRT(PTKEM) * (-2./3.) * PTP
+ELSE
+ ZFLX(:,:,:)= (2./3.) * PTKEM &
+ - XCMFS * PK *( (4./3.) * GZ_W_M_PWM &
+ -(2./3.) * ( GX_U_M_PUM ) )
+ !! & to be tested
+ !! -2.* XCMFB * PLM / SQRT(PTKEM) * (-2./3.) * PTP
+END IF
+!
+ZFLX(:,:,IKE+1)= ZFLX(:,:,IKE)
+!
+ZFLX(:,:,IKB) = (2./3.) * PTKEM(:,:,IKB) &
+ - XCMFS * PK(:,:,IKB) * 2. * ZDW_DZ(:,:,1)
+
+! & to be tested
+! - 2.* XCMFB * PLM(:,:,IKB:IKB) /SQRT(PTKEM(:,:,IKB:IKB)) * &
+! (-2./3.) * PTP(:,:,IKB:IKB)
+!
+! extrapolates this flux under the ground with the surface flux
+ZFLX(:,:,IKB-1) = &
+ PTAU11M(:,:) * ZDIRSINZW(:,:)**2 &
+ + PTAU33M(:,:) * PDIRCOSZW(:,:)**2 &
+ +2. * PCDUEFF(:,:)* PUSLOPEM(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
+ !
+ZFLX(:,:,IKB-1) = 2. * ZFLX(:,:,IKB-1) - ZFLX(:,:,IKB)
+!
+IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+ ! stores <W W>
+ TZFIELD%CMNHNAME = 'W_VAR'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'W_VAR'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_W_VAR'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+END IF
+!
+! Complete the W tendency
+!
+!PRWS(:,:,:)=PRWS(:,:,:) - DZM( PRHODJ*ZFLX/MZF(PDZZ) )
+ZDFDDWDZ(:,:,:) = - XCMFS * PK(:,:,:) * (4./3.)
+ZDFDDWDZ(:,:,:IKB) = 0.
+!
+CALL TRIDIAG_W(PWM,ZFLX,ZDFDDWDZ,PTSTEP,ZMZF_DZZ,PRHODJ,ZWP)
+!
+PRWS = PRWS(:,:,:) + MZM(PRHODJ(:,:,:))*(ZWP(:,:,:)-PWM(:,:,:))/PTSTEP
+!
+!* recomputes flux using guess of W
+!
+GZ_W_M_ZWP = GZ_W_M(ZWP,PDZZ)
+ZFLX(:,:,IKB+1:)=ZFLX(:,:,IKB+1:) &
+ - XCMFS * PK(:,:,IKB+1:) * (4./3.) * (GZ_W_M_ZWP(:,:,IKB+1:) - GZ_W_M_PWM(:,:,IKB+1:))
+!
+IF (KSPLT==1) THEN
+ !Contribution to the dynamic production of TKE:
+! ZWORK(:,:,:) = - ZFLX(:,:,:) * GZ_W_M_PWM
+ ZWORK(:,:,:) = - ZFLX(:,:,:) * GZ_W_M_ZWP
+ !
+ ! evaluate the dynamic production at w(IKB+1) in PDP(IKB)
+ !
+ ZWORK(:,:,IKB) = 0.5* ( -ZFLX(:,:,IKB)*ZDW_DZ(:,:,1) + ZWORK(:,:,IKB+1) )
+ !
+ PDP(:,:,:) = PDP(:,:,:) + ZWORK(:,:,:)
+END IF
+!
+! Storage in the LES configuration
+!
+!
+IF (LLES_CALL .AND. KSPLT==1) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( ZFLX, X_LES_SUBGRID_W2 )
+ CALL LES_MEAN_SUBGRID( -ZWORK, X_LES_RES_ddxa_W_SBG_UaW , .TRUE.)
+ CALL LES_MEAN_SUBGRID( GZ_M_M(PTHLM,PDZZ)*ZFLX, X_LES_RES_ddxa_Thl_SBG_UaW , .TRUE.)
+ CALL LES_MEAN_SUBGRID(ZFLX*MZF(GZ_M_W(1,IKU,1,PTHLM,PDZZ)),X_LES_RES_ddz_Thl_SBG_W2)
+ IF (KRR>=1) THEN
+ CALL LES_MEAN_SUBGRID( GZ_M_M(PRM(:,:,:,1),PDZZ)*ZFLX, &
+ X_LES_RES_ddxa_Rt_SBG_UaW , .TRUE.)
+ CALL LES_MEAN_SUBGRID(ZFLX*MZF(GZ_M_W(1,IKU,1,PRM(:,:,:,1),PDZZ)), &
+ X_LES_RES_ddz_Rt_SBG_W2)
+ END IF
+ DO JSV=1,NSV
+ CALL LES_MEAN_SUBGRID( GZ_M_M(PSVM(:,:,:,JSV),PDZZ)*ZFLX, &
+ X_LES_RES_ddxa_Sv_SBG_UaW(:,:,:,JSV) , .TRUE.)
+ CALL LES_MEAN_SUBGRID(ZFLX*MZF(GZ_M_W(1,IKU,1,PSVM(:,:,:,JSV),PDZZ)), &
+ X_LES_RES_ddz_Sv_SBG_W2(:,:,:,JSV))
+ END DO
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+!
+CALL CLEANLIST_ll(TZFIELDS_ll)
+!
+!
+END SUBROUTINE TURB_HOR_DYN_CORR
+END MODULE MODE_TURB_HOR_DYN_CORR
diff --git a/src/common/turb/mode_turb_hor_splt.F90 b/src/common/turb/mode_turb_hor_splt.F90
new file mode 100644
index 0000000000000000000000000000000000000000..456abc075e74228bce42073ded8d1a332d760d44
--- /dev/null
+++ b/src/common/turb/mode_turb_hor_splt.F90
@@ -0,0 +1,535 @@
+!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_HOR_SPLT
+IMPLICIT NONE
+CONTAINS
+ SUBROUTINE TURB_HOR_SPLT(KSPLIT, KRR, KRRL, KRRI, PTSTEP, &
+ HLBCX,HLBCY,OTURB_FLX,OSUBG_COND, &
+ TPFILE, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, &
+ PDIRCOSXW,PDIRCOSYW,PDIRCOSZW, &
+ PCOSSLOPE,PSINSLOPE, &
+ PRHODJ,PTHVREF, &
+ PSFTHM,PSFRM,PSFSVM, &
+ PCDUEFF,PTAU11M,PTAU12M,PTAU22M,PTAU33M, &
+ PUM,PVM,PWM,PUSLOPEM,PVSLOPEM,PTHLM,PRM,PSVM, &
+ PTKEM,PLM,PLEPS, &
+ PLOCPEXNM,PATHETA,PAMOIST,PSRCM,PFRAC_ICE, &
+ PDP,PTP,PSIGS, &
+ PTRH, &
+ PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS )
+! ################################################################
+!
+!
+!!**** *TURB_HOR* -routine to compute the source terms in the meso-NH
+!! model equations due to the non-vertical turbulent fluxes.
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to compute the non-vertical
+! turbulent fluxes of the evolutive variables and give back the
+! source terms to the main program.
+!
+!!** METHOD
+!! ------
+!! Complementary 3D calculations when running at high resolution;
+!! The non-vertical turbulent fluxes are computed explicitly. The
+!! contributions are cumulated in PRvarS and in DP and TP of TKE
+!
+! d(rho*T) = -d(rho*u'T'/dxx) -d(-rho*u'T'*dzx/dxx/dzz)
+! / dt / dx /dz
+!!
+!!
+!! Near the bottom of the model, uncentred evaluation of vertical
+!! gradients are required because no field values are available under
+!! the level where the gradient must be evaluated. In this case, the
+!! gradient is computed with a second order accurate uncentred scheme
+!! according to:
+!!
+!! D FF dzz3 (dzz3+dzz4)
+!! ---- = - ----------------- FF(4) + ----------------- FF(3)
+!! D z (dzz3+dzz4) dzz4 dzz3 dzz4
+!!
+!! dzz4 + 2 dzz3
+!! - ----------------- FF(2)
+!! (dzz3+dzz4) dzz3
+!!
+!! where the values are taken from:
+!!
+!! ----- FF(5)
+!! |
+!! | dzz5
+!! |
+!! ----- FF(4)
+!! |
+!! | dzz4
+!! |
+!! ----- FF(3)
+!! |
+!! | dzz3
+!! |
+!! ----- FF(2) , (D FF / DZ)
+!! | dzz2 * 0.5
+!! ----- ground
+!!
+!!
+!!
+!! EXTERNAL
+!! --------
+!! GX_M_U, GY_M_V
+!! GX_M_M, GY_M_M, GZ_M_M
+!! GY_U_UV,GX_V_UV
+!! GX_U_M, GY_V_M, GZ_W_M
+!! GX_W_UW,GY_W_UW
+!! : Cartesian vertical gradient operators
+!!
+!!
+!! MXM,MXF,MYM,MYF,MZM,MZF
+!! : Shuman functions (mean operators)
+!! DXM,DXF.DYM,DYF,DZM,DZF
+!! : Shuman functions (difference operators)
+!!
+!!
+!! 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 : number of vertical external points
+!!
+!! Module MODD_CONF
+!!
+!! CPROGRAM
+!!
+!!
+!! REFERENCE
+!! ---------
+!! Book 2 of documentation (routine TURB_HOR)
+!! Book 1 of documentation (Chapter: Turbulence)
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original Aug 29, 1994
+!! Modifications: Feb 14, 1995 (J.Cuxart and J.Stein)
+!! Doctorization and Optimization
+!! March 21, 1995 (J.M. Carriere)
+!! Introduction of cloud water
+!! June 14, 1995 (J. Stein)
+!! rm the ZVTPV computation + bug in the all
+!! or nothing condens. case
+!! June 28, 1995 (J.Cuxart) Add the LES tools
+!! Sept 19, 1995 (J. Stein) change the surface flux
+!! computations
+!! Nov 13, 1995 (J. Stein) include the tangential fluxes
+!! bug in <u'w'> at the surface
+!! Nov 27, 1997 (V. Saravane) spliting of the routine
+!! Nov 27, 1997 (V. Masson) clearing of the routine
+!! Mar 07, 2001 (V. Masson and J. Stein) time splitting
+!! + major bugs correction for slopes
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! Feb 20, 2003 (JP Pinty) Add PFRAC_ICE
+!! Oct.2009 (C.Lac) Introduction of different PTSTEP according to the
+!! advection schemes
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CONF
+USE MODD_CST
+USE MODD_CTURB
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_PARAMETERS
+!
+!
+USE MODI_SHUMAN
+USE MODE_TURB_HOR
+USE MODE_TURB_HOR_TKE
+!
+USE MODE_ll
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declaration of arguments
+!
+!
+INTEGER, INTENT(IN) :: KSPLIT ! number of time splitting
+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.
+REAL, INTENT(IN) :: PTSTEP ! timestep
+CHARACTER (LEN=*), DIMENSION(:), INTENT(IN) :: HLBCX,HLBCY
+LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
+ ! turbulent fluxes in the syncronous FM-file
+LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for sub-grid
+! condensation
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! vertical grid
+REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSXW, PDIRCOSYW, PDIRCOSZW
+! Director Cosinus along x, y and z directions at surface w-point
+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 ! density * grid volume
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! ref. state VPT
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHM,PSFRM
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSFSVM ! surface fluxes
+!
+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) :: PTAU22M ! <vv> in the same axes
+REAL, DIMENSION(:,:), INTENT(IN) :: PTAU33M ! <ww> in the same axes
+!
+! Variables at t-1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUM,PVM,PWM,PTHLM
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! mixing ratios at t-1,
+ ! where PRM(:,:,:,1) = conservative mixing ratio
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! scalar var. at t-1
+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- dt
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLEPS ! dissipative length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLOCPEXNM ! Lv(T)/Cp/Exner 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(INOUT) :: PRUS, PRVS, PRWS, PRTHLS
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRSVS,PRRS ! var. at t+1 -split-
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDP,PTP ! TKE production terms
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PTRH
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSIGS
+ ! IN: Vertical part of Sigma_s at t
+ ! OUT: Total Sigma_s at t
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+REAL,ALLOCATABLE,DIMENSION(:,:,:) :: ZK ! Turbulent diffusion doef.
+ ! ZK = PLM * SQRT(PTKEM)
+REAL,ALLOCATABLE,DIMENSION(:,:,:) :: ZINV_PDXX ! 1./PDXX
+REAL,ALLOCATABLE,DIMENSION(:,:,:) :: ZINV_PDYY ! 1./PDYY
+REAL,ALLOCATABLE,DIMENSION(:,:,:) :: ZINV_PDZZ ! 1./PDZZ
+REAL,ALLOCATABLE,DIMENSION(:,:,:) :: ZMZM_PRHODJ ! MZM(PRHODJ)
+!
+INTEGER :: JSPLT ! current split
+!
+INTEGER :: IKB, IKE, IIB, IIE, IJB, IJE
+INTEGER :: JRR, JSV
+!
+INTEGER :: ISV
+INTEGER :: IINFO_ll
+!
+REAL,ALLOCATABLE,DIMENSION(:,:,:) :: ZUM, ZVM, ZWM, ZTHLM, ZTKEM
+REAL,ALLOCATABLE,DIMENSION(:,:,:,:) :: ZRM, ZSVM
+REAL,ALLOCATABLE,DIMENSION(:,:,:) :: ZRUS, ZRVS, ZRWS, ZRTHLS
+REAL,ALLOCATABLE,DIMENSION(:,:,:,:) :: ZRRS, ZRSVS
+!
+!
+TYPE(LIST_ll), POINTER, SAVE :: TZFIELDS_ll
+!
+! ---------------------------------------------------------------------------
+!
+!* 1. PRELIMINARY COMPUTATIONS
+! ------------------------
+!
+IKB = 1.+JPVEXT
+IKE = SIZE(PUM,3) - JPVEXT
+CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
+ISV=SIZE(PSVM,4)
+!
+ALLOCATE(ZK(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)))
+ALLOCATE(ZINV_PDXX(SIZE(PDXX,1),SIZE(PDXX,2),SIZE(PDXX,3)))
+ALLOCATE(ZINV_PDYY(SIZE(PDYY,1),SIZE(PDYY,2),SIZE(PDYY,3)))
+ALLOCATE(ZINV_PDZZ(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)))
+ALLOCATE(ZMZM_PRHODJ(SIZE(PRHODJ,1),SIZE(PRHODJ,2),SIZE(PRHODJ,3)))
+!
+ZINV_PDXX = 1./PDXX
+ZINV_PDYY = 1./PDYY
+ZINV_PDZZ = 1./PDZZ
+ZMZM_PRHODJ = MZM(PRHODJ)
+!
+ZK(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:))
+!
+NULLIFY(TZFIELDS_ll)
+!
+!--------------------------------------------------------------------
+!
+!* 2. SPLIT PROCESS LOOP
+! ------------------
+!
+IF (KSPLIT>1 .AND. CPROGRAM=='MESONH') THEN
+!
+!* 2.1 allocations
+! -----------
+!
+ ALLOCATE(ZUM(SIZE(PUM,1),SIZE(PUM,2),SIZE(PUM,3)))
+ ALLOCATE(ZVM(SIZE(PVM,1),SIZE(PVM,2),SIZE(PVM,3)))
+ ALLOCATE(ZWM(SIZE(PWM,1),SIZE(PWM,2),SIZE(PWM,3)))
+ ALLOCATE(ZSVM(SIZE(PSVM,1),SIZE(PSVM,2),SIZE(PSVM,3),SIZE(PSVM,4)))
+ ALLOCATE(ZTHLM(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)))
+ ALLOCATE(ZTKEM(SIZE(PTKEM,1),SIZE(PTKEM,2),SIZE(PTKEM,3)))
+ ALLOCATE(ZRM(SIZE(PRM,1),SIZE(PRM,2),SIZE(PRM,3),SIZE(PRM,4)))
+ ALLOCATE(ZRUS(SIZE(PRUS,1),SIZE(PRUS,2),SIZE(PRUS,3)))
+ ALLOCATE(ZRVS(SIZE(PRVS,1),SIZE(PRVS,2),SIZE(PRVS,3)))
+ ALLOCATE(ZRWS(SIZE(PRWS,1),SIZE(PRWS,2),SIZE(PRWS,3)))
+ ALLOCATE(ZRSVS(SIZE(PRSVS,1),SIZE(PRSVS,2),SIZE(PRSVS,3),SIZE(PRSVS,4)))
+ ALLOCATE(ZRTHLS(SIZE(PRTHLS,1),SIZE(PRTHLS,2),SIZE(PRTHLS,3)))
+ ALLOCATE(ZRRS(SIZE(PRRS,1),SIZE(PRRS,2),SIZE(PRRS,3),SIZE(PRRS,4)))
+!
+!
+!* 2.2 list for parallel exchanges
+! ---------------------------
+!
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZUM, 'TURB_HOR_SPLT::ZUM' )
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZVM, 'TURB_HOR_SPLT::ZVM' )
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZWM, 'TURB_HOR_SPLT::ZWM' )
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZTHLM, 'TURB_HOR_SPLT::ZTHLM' )
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZTKEM, 'TURB_HOR_SPLT::ZTKEM' )
+ CALL ADD4DFIELD_ll( TZFIELDS_ll, ZSVM(:,:,:,1:ISV), 'TURB_HOR_SPLT::ZSVM(:,:,:,1:ISV)' )
+ CALL ADD4DFIELD_ll( TZFIELDS_ll, ZRM(:,:,:,1:KRR), 'TURB_HOR_SPLT::ZRM(:,:,:,1:KRR)' )
+!
+!
+!* 2.3 initializations
+! ---------------
+!
+!
+ ZUM=PUM
+ ZVM=PVM
+ ZWM=PWM
+ IF (ISV>0) ZSVM=PSVM
+ ZTHLM=PTHLM
+ ZTKEM=PTKEM
+ IF (KRR>0) ZRM=PRM
+ !
+ ZRUS=PRUS*KSPLIT
+ ZRVS=PRVS*KSPLIT
+ ZRWS=PRWS*KSPLIT
+ IF (ISV>0) ZRSVS=PRSVS*KSPLIT
+ ZRTHLS=PRTHLS*KSPLIT
+ IF (KRR>0) ZRRS=PRRS*KSPLIT
+
+!
+!* 2.4 split process
+! -------------
+!
+ DO JSPLT=1,KSPLIT
+!
+! compute the turbulent tendencies for the small time step
+ CALL TURB_HOR(JSPLT, KRR, KRRL, KRRI, PTSTEP, &
+ OTURB_FLX,OSUBG_COND, &
+ TPFILE, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, &
+ PDIRCOSXW,PDIRCOSYW,PDIRCOSZW, &
+ PCOSSLOPE,PSINSLOPE, &
+ ZINV_PDXX, ZINV_PDYY, ZINV_PDZZ, ZMZM_PRHODJ, &
+ ZK, &
+ PRHODJ,PTHVREF, &
+ PSFTHM,PSFRM,PSFSVM, &
+ PCDUEFF,PTAU11M,PTAU12M,PTAU22M,PTAU33M, &
+ ZUM,ZVM,ZWM,PUSLOPEM,PVSLOPEM,ZTHLM,ZRM,ZSVM, &
+ PTKEM,PLM,PLEPS, &
+ PLOCPEXNM,PATHETA,PAMOIST,PSRCM,PFRAC_ICE, &
+ PDP,PTP,PSIGS, &
+ ZRUS,ZRVS,ZRWS,ZRTHLS,ZRRS,ZRSVS )
+!
+! horizontal transport of Tke
+!
+ CALL TURB_HOR_TKE(JSPLT, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ ZINV_PDXX, ZINV_PDYY, ZINV_PDZZ, ZMZM_PRHODJ, &
+ ZK, PRHODJ, ZTKEM, &
+ PTRH )
+!
+!
+! split temporal advance
+
+ ZUM=PUM+(ZRUS/KSPLIT-PRUS)/MXM(PRHODJ)*PTSTEP
+ ZVM=PVM+(ZRVS/KSPLIT-PRVS)/MYM(PRHODJ)*PTSTEP
+ ZWM=PWM+(ZRWS/KSPLIT-PRWS)/ZMZM_PRHODJ*PTSTEP
+ DO JSV=1,ISV
+ ZSVM(:,:,:,JSV)=PSVM(:,:,:,JSV)+ &
+ (ZRSVS(:,:,:,JSV)/KSPLIT-PRSVS(:,:,:,JSV))/PRHODJ*PTSTEP
+ END DO
+ ZTHLM=PTHLM+(ZRTHLS/KSPLIT-PRTHLS)/PRHODJ*PTSTEP
+ ZTKEM=ZTKEM+PTRH*PTSTEP/KSPLIT
+ DO JRR=1,KRR
+ ZRM(:,:,:,JRR)=PRM(:,:,:,JRR)+ &
+ (ZRRS(:,:,:,JRR)/KSPLIT-PRRS(:,:,:,JRR))/PRHODJ*PTSTEP
+ END DO
+!
+! reinforce boundary conditions
+!
+ IF (JSPLT<KSPLIT-NHALO+1) CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ !
+ IF ( HLBCX(1) /= "CYCL" .AND. LWEST_ll()) THEN
+ ZUM(IIB ,:,:)=PUM(IIB ,:,:)
+ ZVM(IIB-1,:,:)=PVM(IIB-1,:,:)
+ ZWM(IIB-1,:,:)=PWM(IIB-1,:,:)
+ ZTHLM(IIB-1,:,:)=PTHLM(IIB-1,:,:)
+ ZTKEM(IIB-1,:,:)=PTKEM(IIB-1,:,:)
+ IF (ISV>0) ZSVM(IIB-1,:,:,:)=PSVM(IIB-1,:,:,:)
+ IF (KRR>0) ZRM (IIB-1,:,:,:)=PRM (IIB-1,:,:,:)
+ ENDIF
+ !
+ IF ( HLBCX(2) /= "CYCL" .AND. LEAST_ll()) THEN
+ ZUM(IIE+1,:,:)=PUM(IIE+1,:,:)
+ ZVM(IIE+1,:,:)=PVM(IIE+1,:,:)
+ ZWM(IIE+1,:,:)=PWM(IIE+1,:,:)
+ ZTHLM(IIE+1,:,:)=PTHLM(IIE+1,:,:)
+ ZTKEM(IIE+1,:,:)=PTKEM(IIE+1,:,:)
+ IF (ISV>0) ZSVM(IIE+1,:,:,:)=PSVM(IIE+1,:,:,:)
+ IF (KRR>0) ZRM (IIE+1,:,:,:)=PRM(IIE+1,:,:,:)
+ ENDIF
+ !
+ IF ( HLBCY(1) /= "CYCL" .AND. LSOUTH_ll()) THEN
+ ZUM(:,IJB-1,:)=PUM(:,IJB-1,:)
+ ZVM(:,IJB ,:)=PVM(:,IJB ,:)
+ ZWM(:,IJB-1,:)=PWM(:,IJB-1,:)
+ ZTHLM(:,IJB-1,:)=PTHLM(:,IJB-1,:)
+ ZTKEM(:,IJB-1,:)=PTKEM(:,IJB-1,:)
+ IF (ISV>0) ZSVM(:,IJB-1,:,:)=PSVM(:,IJB-1,:,:)
+ IF (KRR>0) ZRM (:,IJB-1,:,:)=PRM (:,IJB-1,:,:)
+ ENDIF
+ !
+ IF ( HLBCY(2) /= "CYCL" .AND. LNORTH_ll()) THEN
+ ZUM(:,IJE+1,:)=PUM(:,IJE+1,:)
+ ZVM(:,IJE+1,:)=PVM(:,IJE+1,:)
+ ZWM(:,IJE+1,:)=PWM(:,IJE+1,:)
+ ZTHLM(:,IJE+1,:)=PTHLM(:,IJE+1,:)
+ ZTKEM(:,IJE+1,:)=PTKEM(:,IJE+1,:)
+ IF (ISV>0) ZSVM(:,IJE+1,:,:)=PSVM(:,IJE+1,:,:)
+ IF (KRR>0) ZRM (:,IJE+1,:,:)=PRM(:,IJE+1,:,:)
+ ENDIF
+ !
+ ZUM(:,:,IKB-1)=ZUM(:,:,IKB)
+ ZVM(:,:,IKB-1)=ZVM(:,:,IKB)
+ ZWM(:,:,IKB-1)=ZWM(:,:,IKB)
+ ZTHLM(:,:,IKB-1)=ZTHLM(:,:,IKB)
+ ZTKEM(:,:,IKB-1)=ZTKEM(:,:,IKB)
+ IF (ISV>0) ZSVM(:,:,IKB-1,:)=ZSVM(:,:,IKB,:)
+ IF (KRR>0) ZRM (:,:,IKB-1,:)=ZRM (:,:,IKB,:)
+ !
+ ZUM(:,:,IKE+1)=ZUM(:,:,IKE)
+ ZVM(:,:,IKE+1)=ZVM(:,:,IKE)
+ ZWM(:,:,IKE+1)=ZWM(:,:,IKE)
+ ZTHLM(:,:,IKE+1)=ZTHLM(:,:,IKE)
+ ZTKEM(:,:,IKE+1)=ZTKEM(:,:,IKE)
+ IF (ISV>0) ZSVM(:,:,IKE+1,:)=ZSVM(:,:,IKE,:)
+ IF (KRR>0) ZRM (:,:,IKE+1,:)=ZRM (:,:,IKE,:)
+ !
+ END DO
+!
+!* 2.5 update the complete tendencies
+! ------------------------------
+!
+ PRUS=ZRUS/KSPLIT
+ PRVS=ZRVS/KSPLIT
+ PRWS=ZRWS/KSPLIT
+ IF (ISV>0) PRSVS=ZRSVS/KSPLIT
+ PRTHLS=ZRTHLS/KSPLIT
+ IF (KRR>0) PRRS=ZRRS/KSPLIT
+ PTRH=(ZTKEM-PTKEM)/PTSTEP
+!
+!* 2.6 deallocations
+! -------------
+!
+ DEALLOCATE(ZUM)
+ DEALLOCATE(ZVM)
+ DEALLOCATE(ZWM)
+ DEALLOCATE(ZSVM)
+ DEALLOCATE(ZTHLM)
+ DEALLOCATE(ZTKEM)
+ DEALLOCATE(ZRM)
+ DEALLOCATE(ZRUS)
+ DEALLOCATE(ZRVS)
+ DEALLOCATE(ZRWS)
+ DEALLOCATE(ZRSVS)
+ DEALLOCATE(ZRTHLS)
+ DEALLOCATE(ZRRS)
+ !
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+!
+!-------------------------------------------------------------------
+!
+!* 4. NO SPLIT PROCESS CASE
+! ---------------------
+!
+ELSE
+!
+ CALL TURB_HOR(1, KRR, KRRL, KRRI, PTSTEP, &
+ OTURB_FLX,OSUBG_COND, &
+ TPFILE, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, &
+ PDIRCOSXW,PDIRCOSYW,PDIRCOSZW, &
+ PCOSSLOPE,PSINSLOPE, &
+ ZINV_PDXX, ZINV_PDYY, ZINV_PDZZ, ZMZM_PRHODJ, &
+ ZK, &
+ PRHODJ,PTHVREF, &
+ PSFTHM,PSFRM,PSFSVM, &
+ PCDUEFF,PTAU11M,PTAU12M,PTAU22M,PTAU33M, &
+ PUM,PVM,PWM,PUSLOPEM,PVSLOPEM,PTHLM,PRM,PSVM, &
+ PTKEM,PLM,PLEPS, &
+ PLOCPEXNM,PATHETA,PAMOIST,PSRCM,PFRAC_ICE, &
+ PDP,PTP,PSIGS, &
+ PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS )
+
+! horizontal transport of Tke
+!
+
+ CALL TURB_HOR_TKE(1, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ ZINV_PDXX, ZINV_PDYY, ZINV_PDZZ, ZMZM_PRHODJ, &
+ ZK, PRHODJ, PTKEM, &
+ PTRH )
+!
+END IF
+!--------------------------------------------------------------------
+!
+DEALLOCATE(ZK)
+DEALLOCATE(ZINV_PDXX)
+DEALLOCATE(ZINV_PDYY)
+DEALLOCATE(ZINV_PDZZ)
+DEALLOCATE(ZMZM_PRHODJ)
+!
+END SUBROUTINE TURB_HOR_SPLT
+END MODULE MODE_TURB_HOR_SPLT
diff --git a/src/common/turb/mode_turb_hor_sv_corr.F90 b/src/common/turb/mode_turb_hor_sv_corr.F90
new file mode 100644
index 0000000000000000000000000000000000000000..ef090e0772dea9b5bf82be69c2d855dd9f2169c2
--- /dev/null
+++ b/src/common/turb/mode_turb_hor_sv_corr.F90
@@ -0,0 +1,185 @@
+!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_HOR_SV_CORR
+IMPLICIT NONE
+CONTAINS
+ SUBROUTINE TURB_HOR_SV_CORR(KRR,KRRL,KRRI, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ PLM,PLEPS,PTKEM,PTHVREF, &
+ PTHLM,PRM, &
+ PLOCPEXNM,PATHETA,PAMOIST,PSRCM, &
+ PWM,PSVM )
+! ################################################################
+!
+!
+!!**** *TURB_HOT_SV_CORR* computes subgrid Sv2 and SvThv terms
+!!
+!! PURPOSE
+!! -------
+!!
+!! see TURB_HOR
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!!
+!! V. Masson * Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 06/11/02
+!! JP Pinty Feb 20, 2003 Add PFRAC_ICE
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+USE MODD_CONF
+USE MODD_CTURB
+USE MODD_PARAMETERS
+USE MODD_NSV, ONLY : NSV,NSV_LGBEG,NSV_LGEND
+USE MODD_LES
+USE MODD_BLOWSNOW
+!
+USE MODI_GRADIENT_M
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_SHUMAN
+USE MODI_LES_MEAN_SUBGRID
+USE MODE_EMOIST, ONLY: EMOIST
+USE MODE_ETHETA, ONLY: ETHETA
+!
+USE MODI_SECOND_MNH
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declaration of arguments
+!
+!
+!
+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) :: PDXX, PDYY, PDZZ, PDZX, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! mixing length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLEPS ! dissipative length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! tke
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! reference Thv
+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) :: 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) :: PWM ! w at t-1
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! scalar var. at t-1
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+REAL, DIMENSION(SIZE(PSVM,1),SIZE(PSVM,2),SIZE(PSVM,3)) &
+ :: ZFLX, ZA
+!
+INTEGER :: JSV ! loop counter
+INTEGER :: IKU
+!
+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 :: ZCSV !constant for the scalar flux
+! ---------------------------------------------------------------------------
+!
+IKU=SIZE(PTKEM,3)
+CALL SECOND_MNH(ZTIME1)
+!
+IF(LBLOWSNOW) THEN
+! See Vionnet (PhD, 2012) for a complete discussion around the value of the Schmidt number for blowing snow variables
+ ZCSV= XCHF/XRSNOW
+ELSE
+ ZCSV= XCHF
+ENDIF
+!
+DO JSV=1,NSV
+!
+ IF (LNOMIXLG .AND. JSV >= NSV_LGBEG .AND. JSV<= NSV_LGEND) CYCLE
+ !
+ ! variance Sv2
+ !
+ IF (LLES_CALL) THEN
+ IF (.NOT. L2D) THEN
+ ZFLX(:,:,:) = ZCSV / ZCSVD * PLM(:,:,:) * PLEPS(:,:,:) * &
+ ( GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)**2 &
+ + GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)**2 )
+ ELSE
+ ZFLX(:,:,:) = ZCSV / ZCSVD * PLM(:,:,:) * PLEPS(:,:,:) * &
+ GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)**2
+ END IF
+ CALL LES_MEAN_SUBGRID( -2.*ZCSVD*SQRT(PTKEM)*ZFLX/PLEPS, &
+ X_LES_SUBGRID_DISS_Sv2(:,:,:,JSV), .TRUE. )
+ CALL LES_MEAN_SUBGRID( MZF(PWM)*ZFLX, X_LES_RES_W_SBG_Sv2(:,:,:,JSV), .TRUE. )
+ END IF
+ !
+ ! covariance SvThv
+ !
+ IF (LLES_CALL) THEN
+ ZA(:,:,:) = ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM)
+ IF (.NOT. L2D) THEN
+ ZFLX(:,:,:)= PLM(:,:,:) * PLEPS(:,:,:) &
+ * ( GX_M_M(PTHLM,PDXX,PDZZ,PDZX) * GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX) &
+ + GY_M_M(PTHLM,PDYY,PDZZ,PDZY) * GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY) &
+ ) * (XCSHF+ZCSV) / (2.*ZCTSVD)
+ ELSE
+ ZFLX(:,:,:)= PLM(:,:,:) * PLEPS(:,:,:) &
+ * GX_M_M(PTHLM,PDXX,PDZZ,PDZX) * GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX) &
+ * (XCSHF+ZCSV) / (2.*ZCTSVD)
+ END IF
+ CALL LES_MEAN_SUBGRID( ZA*ZFLX, X_LES_SUBGRID_SvThv(:,:,:,JSV) , .TRUE.)
+ CALL LES_MEAN_SUBGRID( -XG/PTHVREF/3.*ZA*ZFLX, X_LES_SUBGRID_SvPz(:,:,:,JSV), .TRUE. )
+ !
+ IF (KRR>=1) THEN
+ ZA(:,:,:) = EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM)
+ IF (.NOT. L2D) THEN
+ ZFLX(:,:,:)= PLM(:,:,:) * PLEPS(:,:,:) &
+ * ( GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX) * GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX) &
+ + GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY) * GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY) &
+ ) * (XCHF+ZCSV) / (2.*ZCQSVD)
+ ELSE
+ ZFLX(:,:,:)= PLM(:,:,:) * PLEPS(:,:,:) &
+ * GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX) * GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX) &
+ * (XCHF+ZCSV) / (2.*ZCQSVD)
+ END IF
+ CALL LES_MEAN_SUBGRID( ZA*ZFLX, X_LES_SUBGRID_SvThv(:,:,:,JSV) , .TRUE.)
+ CALL LES_MEAN_SUBGRID( -XG/PTHVREF/3.*ZA*ZFLX, X_LES_SUBGRID_SvPz(:,:,:,JSV), .TRUE. )
+ END IF
+ END IF
+!
+END DO ! end loop JSV
+!
+CALL SECOND_MNH(ZTIME2)
+XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+!
+END SUBROUTINE TURB_HOR_SV_CORR
+END MODULE MODE_TURB_HOR_SV_CORR
+
diff --git a/src/common/turb/mode_turb_hor_sv_flux.F90 b/src/common/turb/mode_turb_hor_sv_flux.F90
new file mode 100644
index 0000000000000000000000000000000000000000..f61a531d910f2c0fb1f97b249a8665a53fa1ca08
--- /dev/null
+++ b/src/common/turb/mode_turb_hor_sv_flux.F90
@@ -0,0 +1,315 @@
+!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_HOR_SV_FLUX
+IMPLICIT NONE
+CONTAINS
+ SUBROUTINE TURB_HOR_SV_FLUX(KSPLT, &
+ OTURB_FLX, &
+ TPFILE, &
+ PK,PINV_PDXX,PINV_PDYY,PINV_PDZZ,PMZM_PRHODJ, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ PDIRCOSXW,PDIRCOSYW, &
+ PRHODJ,PWM, &
+ PSFSVM, &
+ PSVM, &
+ PRSVS )
+! ################################################################
+!
+!
+!!**** *TURB_HOR* -routine to compute the source terms in the meso-NH
+!! model equations due to the non-vertical turbulent fluxes.
+!!
+!! PURPOSE
+!! -------
+!!
+!! see TURB_HOR
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!!
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Aug , 1997 (V. Saravane) spliting of TURB_HOR
+!! Nov 27, 1997 (V. Masson) clearing of the routine
+!! Oct 18, 2000 (V. Masson) LES computations + LFLAT swith
+!! + bug on Y scalar flux
+!! 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
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+USE MODD_CONF
+USE MODD_CTURB
+USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_PARAMETERS
+USE MODD_NSV, ONLY: NSV_LGBEG, NSV_LGEND
+USE MODD_LES
+USE MODD_BLOWSNOW
+!
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
+!
+USE MODI_GRADIENT_M
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_SHUMAN
+USE MODE_COEFJ, ONLY: COEFJ
+USE MODI_LES_MEAN_SUBGRID
+!
+USE MODI_SECOND_MNH
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declaration of arguments
+!
+!
+!
+INTEGER, INTENT(IN) :: KSPLT ! split process index
+LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
+ ! turbulent fluxes in the syncronous FM-file
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PK ! Turbulent diffusion doef.
+ ! PK = PLM * SQRT(PTKEM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDXX ! 1./PDXX
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDYY ! 1./PDYY
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDZZ ! 1./PDZZ
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PMZM_PRHODJ ! MZM(PRHODJ)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSXW, PDIRCOSYW
+! Director Cosinus along x and y directions at surface w-point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density * grid volume
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM ! vertical wind
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSFSVM ! surface fluxes
+!
+!
+! Variables at t-1
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! scalar var. at t-1
+!
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRSVS ! var. at t+1 -split-
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+REAL, DIMENSION(SIZE(PSVM,1),SIZE(PSVM,2),SIZE(PSVM,3)) &
+ :: ZFLXX,ZFLXY
+ ! work arrays
+REAL, DIMENSION(SIZE(PSVM,1),SIZE(PSVM,2),1) :: ZWORK2D
+!
+REAL :: ZCSV !constant for the scalar flux
+
+INTEGER :: IKB,IKE
+ ! Index values for the Beginning and End
+ ! mass points of the domain
+INTEGER :: JSV ! loop counter
+INTEGER :: ISV ! number of scalar var.
+REAL, DIMENSION(SIZE(PDZZ,1),SIZE(PDZZ,2),1+JPVEXT:3+JPVEXT) :: ZCOEFF
+ ! coefficients for the uncentred gradient
+ ! computation near the ground
+!
+INTEGER :: IKU
+TYPE(TFIELDDATA) :: TZFIELD
+REAL :: ZTIME1, ZTIME2
+! ---------------------------------------------------------------------------
+!
+!* 1. PRELIMINARY COMPUTATIONS
+! ------------------------
+!
+IKB = 1+JPVEXT
+IKE = SIZE(PSVM,3)-JPVEXT
+IKU = SIZE(PSVM,3)
+!
+ISV = SIZE(PSVM,4)
+!
+IF(LBLOWSNOW) THEN
+! See Vionnet (PhD, 2012) for a complete discussion around the value of the Schmidt number for blowing snow variables
+ ZCSV= XCHF/XRSNOW
+ELSE
+ ZCSV= XCHF
+ENDIF
+!
+! compute the coefficients for the uncentred gradient computation near the
+! ground
+ZCOEFF(:,:,IKB+2)= - PDZZ(:,:,IKB+1) / &
+ ( (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) * PDZZ(:,:,IKB+2) )
+ZCOEFF(:,:,IKB+1)= (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) / &
+ ( PDZZ(:,:,IKB+1) * PDZZ(:,:,IKB+2) )
+ZCOEFF(:,:,IKB)= - (PDZZ(:,:,IKB+2)+2.*PDZZ(:,:,IKB+1)) / &
+ ( (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) * PDZZ(:,:,IKB+1) )
+!
+!
+!* 15. HORIZONTAL FLUXES OF PASSIVE SCALARS
+! ------------------------------------
+!
+!
+DO JSV=1,ISV
+!
+ IF (LNOMIXLG .AND. JSV >= NSV_LGBEG .AND. JSV<= NSV_LGEND) CYCLE
+!
+! 15.1 <U' SVth'>
+! ----------
+!
+ ! Computes the flux in the X direction
+ ZFLXX(:,:,:) = -ZCSV * MXM(PK) * GX_M_U(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)
+ ZFLXX(:,:,IKE+1) = ZFLXX(:,:,IKE)
+!
+! Compute the flux at the first inner U-point with an uncentred vertical
+! gradient
+ ZFLXX(:,:,IKB:IKB) = -ZCSV * MXM( PK(:,:,IKB:IKB) ) * &
+ ( DXM(PSVM(:,:,IKB:IKB,JSV)) * PINV_PDXX(:,:,IKB:IKB) &
+ -MXM ( ZCOEFF(:,:,IKB+2:IKB+2)*PSVM(:,:,IKB+2:IKB+2,JSV) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PSVM(:,:,IKB+1:IKB+1,JSV) &
+ +ZCOEFF(:,:,IKB :IKB )*PSVM(:,:,IKB :IKB ,JSV) &
+ ) * 0.5 * ( PDZX(:,:,IKB+1:IKB+1)+PDZX(:,:,IKB:IKB) ) &
+ * PINV_PDXX(:,:,IKB:IKB) &
+ )
+! extrapolates the flux under the ground so that the vertical average with
+! the IKB flux gives the ground value
+ ZWORK2D(:,:,1)=PSFSVM(:,:,JSV) * PDIRCOSXW(:,:)
+ ZFLXX(:,:,IKB-1:IKB-1) = 2. * MXM( ZWORK2D(:,:,1:1) ) - ZFLXX(:,:,IKB:IKB)
+ !
+ ! stores <U SVth>
+ IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+ WRITE(TZFIELD%CMNHNAME,'("USV_FLX_",I3.3)') JSV
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ TZFIELD%CUNITS = 'SVUNIT m s-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_'//TRIM(TZFIELD%CMNHNAME)
+ TZFIELD%NGRID = 2
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLXX)
+ END IF
+!
+ IF (LLES_CALL .AND. KSPLT==1) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( MXF(ZFLXX), X_LES_SUBGRID_USv(:,:,:,JSV) )
+ CALL LES_MEAN_SUBGRID( MZF(MXF(GX_W_UW(PWM,PDXX,PDZZ,PDZX)*MZM(ZFLXX))), &
+ X_LES_RES_ddxa_W_SBG_UaSv(:,:,:,JSV) , .TRUE. )
+ CALL LES_MEAN_SUBGRID( GX_M_M(PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)*MXF(ZFLXX), &
+ X_LES_RES_ddxa_Sv_SBG_UaSv(:,:,:,JSV), .TRUE. )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+!
+! 15.2 <V' SVth'>
+! ----------
+!
+ IF (.NOT. L2D) THEN
+!
+! Computes the flux in the Y direction
+ ZFLXY(:,:,:)=-ZCSV * MYM(PK) * GY_M_V(1,IKU,1,PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)
+ ZFLXY(:,:,IKE+1) = ZFLXY(:,:,IKE)
+!
+! Compute the flux at the first inner V-point with an uncentred vertical
+! gradient
+!
+ ZFLXY(:,:,IKB:IKB) = -ZCSV * MYM( PK(:,:,IKB:IKB) ) * &
+ ( DYM(PSVM(:,:,IKB:IKB,JSV)) * PINV_PDYY(:,:,IKB:IKB) &
+ -MYM ( ZCOEFF(:,:,IKB+2:IKB+2)*PSVM(:,:,IKB+2:IKB+2,JSV) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PSVM(:,:,IKB+1:IKB+1,JSV) &
+ +ZCOEFF(:,:,IKB :IKB )*PSVM(:,:,IKB :IKB ,JSV) &
+ ) * 0.5 * ( PDZY(:,:,IKB+1:IKB+1)+PDZY(:,:,IKB:IKB) ) &
+ * PINV_PDYY(:,:,IKB:IKB) &
+ )
+! extrapolates the flux under the ground so that the vertical average with
+! the IKB flux gives the ground value
+ ZWORK2D(:,:,1)=PSFSVM(:,:,JSV) * PDIRCOSYW(:,:)
+ ZFLXY(:,:,IKB-1:IKB-1) = 2. * MYM( ZWORK2D(:,:,1:1) ) - ZFLXY(:,:,IKB:IKB)
+ !
+ ! stores <V SVth>
+ IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+ WRITE(TZFIELD%CMNHNAME,'("VSV_FLX_",I3.3)') JSV
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ TZFIELD%CUNITS = 'SVUNIT m s-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_'//TRIM(TZFIELD%CMNHNAME)
+ TZFIELD%NGRID = 3
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLXY)
+ END IF
+!
+ ELSE
+ ZFLXY=0.
+ END IF
+!
+ IF (LLES_CALL .AND. KSPLT==1) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( MYF(ZFLXY), X_LES_SUBGRID_VSv(:,:,:,JSV) )
+ CALL LES_MEAN_SUBGRID( MZF(MYF(GY_W_VW(PWM,PDYY,PDZZ,PDZY)*MZM(ZFLXY))), &
+ X_LES_RES_ddxa_W_SBG_UaSv(:,:,:,JSV) , .TRUE. )
+ CALL LES_MEAN_SUBGRID( GY_M_M(PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)*MYF(ZFLXY), &
+ X_LES_RES_ddxa_Sv_SBG_UaSv(:,:,:,JSV) , .TRUE. )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+!
+!
+! 15.3 Horizontal source terms
+! -----------------------
+!
+ IF (.NOT. L2D) THEN
+ IF (.NOT. LFLAT) THEN
+ PRSVS(:,:,:,JSV)= PRSVS(:,:,:,JSV) &
+ -DXF( MXM(PRHODJ) * ZFLXX * PINV_PDXX ) &
+ -DYF( MYM(PRHODJ) * ZFLXY * PINV_PDYY ) &
+ +DZF( PMZM_PRHODJ * PINV_PDZZ * &
+ ( MXF( MZM(ZFLXX * PINV_PDXX) * PDZX ) + MYF( MZM(ZFLXY * PINV_PDYY) * PDZY ) ) &
+ )
+ ELSE
+ PRSVS(:,:,:,JSV)= PRSVS(:,:,:,JSV) &
+ -DXF( MXM(PRHODJ) * ZFLXX * PINV_PDXX ) &
+ -DYF( MYM(PRHODJ) * ZFLXY * PINV_PDYY )
+ END IF
+ ELSE
+ IF (.NOT. LFLAT) THEN
+ PRSVS(:,:,:,JSV)= PRSVS(:,:,:,JSV) &
+ -DXF( MXM(PRHODJ) * ZFLXX * PINV_PDXX ) &
+ +DZF( PMZM_PRHODJ * PINV_PDZZ * &
+ ( MXF( MZM(ZFLXX * PINV_PDXX) * PDZX ) ) &
+ )
+ ELSE
+ PRSVS(:,:,:,JSV)= PRSVS(:,:,:,JSV) &
+ -DXF( MXM(PRHODJ) *ZFLXX * PINV_PDXX )
+ END IF
+ END IF
+!
+!
+END DO ! end loop JSV
+!
+!
+END SUBROUTINE TURB_HOR_SV_FLUX
+END MODULE MODE_TURB_HOR_SV_FLUX
diff --git a/src/common/turb/mode_turb_hor_thermo_corr.F90 b/src/common/turb/mode_turb_hor_thermo_corr.F90
new file mode 100644
index 0000000000000000000000000000000000000000..f844fbabd7fa660305f976c8797cbecc87c684ea
--- /dev/null
+++ b/src/common/turb/mode_turb_hor_thermo_corr.F90
@@ -0,0 +1,409 @@
+!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_HOR_THERMO_CORR
+IMPLICIT NONE
+CONTAINS
+ SUBROUTINE TURB_HOR_THERMO_CORR(KRR, KRRL, KRRI, &
+ OTURB_FLX,OSUBG_COND, &
+ TPFILE, &
+ PINV_PDXX,PINV_PDYY, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ PTHVREF, &
+ PWM,PTHLM,PRM, &
+ PTKEM,PLM,PLEPS, &
+ PLOCPEXNM,PATHETA,PAMOIST,PSRCM, &
+ PSIGS )
+! ################################################################
+!
+!
+!!**** *TURB_HOR* -routine to compute the source terms in the meso-NH
+!! model equations due to the non-vertical turbulent fluxes.
+!!
+!! PURPOSE
+!! -------
+!!
+!! see TURB_HOR
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!!
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Aug , 1997 (V. Saravane) spliting of TURB_HOR
+!! Nov 27, 1997 (V. Masson) clearing of the routine
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! Feb 20, 2003 (JP Pinty) Add PFRAC_ICE
+!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
+!! change of YCOMMENT
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+USE MODD_CONF
+USE MODD_CTURB
+USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_PARAMETERS
+USE MODD_LES
+!
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
+!
+USE MODI_GRADIENT_M
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_SHUMAN
+USE MODI_LES_MEAN_SUBGRID
+!
+USE MODE_EMOIST, ONLY: EMOIST
+USE MODE_ETHETA, ONLY: ETHETA
+!
+USE MODI_SECOND_MNH
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declaration of arguments
+!
+!
+!
+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
+LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for sub-grid
+! condensation
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDXX ! 1./PDXX
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDYY ! 1./PDYY
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! ref. state Virtual
+ ! Potential Temperature
+!
+! Variables at t-1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! mixing ratios at t-1,
+ ! where PRM(:,:,:,1) = conservative mixing ratio
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! Turb. Kin. Energy
+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
+!
+!
+!
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSIGS
+ ! IN: Vertical part of Sigma_s at t
+ ! OUT: Total Sigma_s at t
+!
+!* 0.2 declaration of local variables
+!
+REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) &
+ :: ZFLX,ZWORK,ZA
+ ! work arrays
+!
+INTEGER :: IKB,IKE
+ ! Index values for the Beginning and End
+ ! mass points of the domain
+REAL, DIMENSION(SIZE(PDZZ,1),SIZE(PDZZ,2),1+JPVEXT:3+JPVEXT) :: ZCOEFF
+ ! coefficients for the uncentred gradient
+ ! computation near the ground
+REAL :: ZTIME1, ZTIME2
+TYPE(TFIELDDATA) :: TZFIELD
+!
+! ---------------------------------------------------------------------------
+!
+!* 1. PRELIMINARY COMPUTATIONS
+! ------------------------
+!
+IKB = 1+JPVEXT
+IKE = SIZE(PTHLM,3)-JPVEXT
+!
+!
+!
+! compute the coefficients for the uncentred gradient computation near the
+! ground
+ZCOEFF(:,:,IKB+2)= - PDZZ(:,:,IKB+1) / &
+ ( (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) * PDZZ(:,:,IKB+2) )
+ZCOEFF(:,:,IKB+1)= (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) / &
+ ( PDZZ(:,:,IKB+1) * PDZZ(:,:,IKB+2) )
+ZCOEFF(:,:,IKB)= - (PDZZ(:,:,IKB+2)+2.*PDZZ(:,:,IKB+1)) / &
+ ( (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) * PDZZ(:,:,IKB+1) )
+!
+!
+!* 8. TURBULENT CORRELATIONS : <THl THl>, <THl Rnp>, <Rnp Rnp>, Sigma_s
+! -----------------------------------------------------------------
+!
+!
+!
+IF ( ( KRRL > 0 .AND. OSUBG_COND) .OR. ( OTURB_FLX .AND. TPFILE%LOPENED ) &
+ .OR. ( LLES_CALL ) ) THEN
+!
+!* 8.1 <THl THl>
+!
+ ! Computes the horizontal variance <THl THl>
+ IF (.NOT. L2D) THEN
+ ZFLX(:,:,:) = XCTV * PLM(:,:,:) * PLEPS(:,:,:) * &
+ ( GX_M_M(PTHLM,PDXX,PDZZ,PDZX)**2 + GY_M_M(PTHLM,PDYY,PDZZ,PDZY)**2 )
+ ELSE
+ ZFLX(:,:,:) = XCTV * PLM(:,:,:) * PLEPS(:,:,:) * &
+ GX_M_M(PTHLM,PDXX,PDZZ,PDZX)**2
+ END IF
+!
+! Compute the flux at the first inner U-point with an uncentred vertical
+! gradient
+!
+ ZFLX(:,:,IKB:IKB) = XCTV * PLM(:,:,IKB:IKB) &
+ * PLEPS(:,:,IKB:IKB) * ( &
+ ( MXF(DXM(PTHLM(:,:,IKB:IKB)) * PINV_PDXX(:,:,IKB:IKB)) &
+ - ( ZCOEFF(:,:,IKB+2:IKB+2)*PTHLM(:,:,IKB+2:IKB+2) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PTHLM(:,:,IKB+1:IKB+1) &
+ +ZCOEFF(:,:,IKB :IKB )*PTHLM(:,:,IKB :IKB ) &
+ ) * 0.5 * ( PDZX(:,:,IKB+1:IKB+1)+PDZX(:,:,IKB:IKB) ) &
+ / MXF(PDXX(:,:,IKB:IKB)) &
+ ) ** 2 + &
+ ( MYF(DYM(PTHLM(:,:,IKB:IKB)) * PINV_PDYY(:,:,IKB:IKB)) &
+ - ( ZCOEFF(:,:,IKB+2:IKB+2)*PTHLM(:,:,IKB+2:IKB+2) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PTHLM(:,:,IKB+1:IKB+1) &
+ +ZCOEFF(:,:,IKB :IKB )*PTHLM(:,:,IKB :IKB ) &
+ ) * 0.5 * ( PDZY(:,:,IKB+1:IKB+1)+PDZY(:,:,IKB:IKB) ) &
+ / MYF(PDYY(:,:,IKB:IKB)) &
+ ) ** 2 )
+ !
+ ZFLX(:,:,IKB-1) = ZFLX(:,:,IKB)
+ !
+ IF ( KRRL > 0 ) THEN
+ ZWORK(:,:,:) = ZFLX(:,:,:) * PATHETA(:,:,:) * PATHETA(:,:,:)
+ END IF
+ !
+ ! stores <THl THl>
+ IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+ TZFIELD%CMNHNAME = 'THL_HVAR'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'THL_HVAR'
+ TZFIELD%CUNITS = 'K2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_THL_HVAR'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+ END IF
+!
+! Storage in the LES configuration (addition to TURB_VER computation)
+!
+ IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( ZFLX, X_LES_SUBGRID_Thl2, .TRUE. )
+ CALL LES_MEAN_SUBGRID( MZF(PWM)*ZFLX, X_LES_RES_W_SBG_Thl2, .TRUE. )
+ CALL LES_MEAN_SUBGRID( -2.*XCTD*SQRT(PTKEM)*ZFLX/PLEPS ,X_LES_SUBGRID_DISS_Thl2, .TRUE. )
+ ZA(:,:,:) = ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM)
+ CALL LES_MEAN_SUBGRID( ZA*ZFLX, X_LES_SUBGRID_ThlThv, .TRUE. )
+ CALL LES_MEAN_SUBGRID( -XG/PTHVREF/3.*ZA*ZFLX, X_LES_SUBGRID_ThlPz, .TRUE. )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+!
+ IF ( KRR /= 0 ) THEN
+!
+!* 8.3 <THl Rnp>
+!
+ ! Computes the horizontal correlation <THl Rnp>
+ IF (.NOT. L2D) THEN
+ ZFLX(:,:,:)= &
+ PLM(:,:,:) * PLEPS(:,:,:) * &
+ (GX_M_M(PTHLM,PDXX,PDZZ,PDZX) * GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX) &
+ + GY_M_M(PTHLM,PDYY,PDZZ,PDZY) * GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY) &
+ ) * (XCHT1+XCHT2)
+ ELSE
+ ZFLX(:,:,:)= &
+ PLM(:,:,:) * PLEPS(:,:,:) * &
+ (GX_M_M(PTHLM,PDXX,PDZZ,PDZX) * GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX) &
+ ) * (XCHT1+XCHT2)
+
+ END IF
+!
+! Compute the flux at the first inner U-point with an uncentred vertical
+! gradient
+ ZFLX(:,:,IKB:IKB) = (XCHT1+XCHT2) * PLM(:,:,IKB:IKB) &
+ * PLEPS(:,:,IKB:IKB) * ( &
+ ( MXF(DXM(PTHLM(:,:,IKB:IKB)) * PINV_PDXX(:,:,IKB:IKB)) &
+ - ( ZCOEFF(:,:,IKB+2:IKB+2)*PTHLM(:,:,IKB+2:IKB+2) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PTHLM(:,:,IKB+1:IKB+1) &
+ +ZCOEFF(:,:,IKB :IKB )*PTHLM(:,:,IKB :IKB ) &
+ ) * 0.5 * ( PDZX(:,:,IKB+1:IKB+1)+PDZX(:,:,IKB:IKB) ) &
+ / MXF(PDXX(:,:,IKB:IKB)) &
+ ) * &
+ ( MXF(DXM(PRM(:,:,IKB:IKB,1)) * PINV_PDXX(:,:,IKB:IKB)) &
+ - ( ZCOEFF(:,:,IKB+2:IKB+2)*PRM(:,:,IKB+2:IKB+2,1) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PRM(:,:,IKB+1:IKB+1,1) &
+ +ZCOEFF(:,:,IKB :IKB )*PRM(:,:,IKB :IKB ,1) &
+ ) * 0.5 * ( PDZX(:,:,IKB+1:IKB+1)+PDZX(:,:,IKB:IKB) ) &
+ / MXF(PDXX(:,:,IKB:IKB)) &
+ ) + &
+ ( MYF(DYM(PTHLM(:,:,IKB:IKB)) * PINV_PDYY(:,:,IKB:IKB)) &
+ - ( ZCOEFF(:,:,IKB+2:IKB+2)*PTHLM(:,:,IKB+2:IKB+2) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PTHLM(:,:,IKB+1:IKB+1) &
+ +ZCOEFF(:,:,IKB :IKB )*PTHLM(:,:,IKB :IKB ) &
+ ) * 0.5 * ( PDZY(:,:,IKB+1:IKB+1)+PDZY(:,:,IKB:IKB) ) &
+ / MYF(PDYY(:,:,IKB:IKB)) &
+ ) * &
+ ( MYF(DYM(PRM(:,:,IKB:IKB,1)) * PINV_PDYY(:,:,IKB:IKB)) &
+ - ( ZCOEFF(:,:,IKB+2:IKB+2)*PRM(:,:,IKB+2:IKB+2,1) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PRM(:,:,IKB+1:IKB+1,1) &
+ +ZCOEFF(:,:,IKB :IKB )*PRM(:,:,IKB :IKB ,1) &
+ ) * 0.5 * ( PDZY(:,:,IKB+1:IKB+1)+PDZY(:,:,IKB:IKB) ) &
+ / MYF(PDYY(:,:,IKB:IKB)) &
+ ) )
+ !
+ ZFLX(:,:,IKB-1) = ZFLX(:,:,IKB)
+ !
+ IF ( KRRL > 0 ) THEN
+ ZWORK(:,:,:) = ZWORK(:,:,:) + &
+ 2. * PATHETA(:,:,:) * PAMOIST(:,:,:) * ZFLX(:,:,:)
+ END IF
+ !
+ ! stores <THl Rnp>
+ IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+ TZFIELD%CMNHNAME = 'THLR_HCOR'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'THLR_HCOR'
+ TZFIELD%CUNITS = 'K kg kg-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_THLR_HCOR'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+ END IF
+!
+! Storage in the LES configuration (addition to TURB_VER computation)
+!
+ IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( ZFLX, X_LES_SUBGRID_ThlRt, .TRUE. )
+ CALL LES_MEAN_SUBGRID( MZF(PWM)*ZFLX, X_LES_RES_W_SBG_ThlRt, .TRUE. )
+ CALL LES_MEAN_SUBGRID( -XCTD*SQRT(PTKEM)*ZFLX/PLEPS ,X_LES_SUBGRID_DISS_ThlRt, .TRUE. )
+ CALL LES_MEAN_SUBGRID( ZA*ZFLX, X_LES_SUBGRID_RtThv, .TRUE. )
+ CALL LES_MEAN_SUBGRID( -XG/PTHVREF/3.*ZA*ZFLX, X_LES_SUBGRID_RtPz,.TRUE.)
+ ZA(:,:,:) = EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM)
+ CALL LES_MEAN_SUBGRID( ZA*ZFLX, X_LES_SUBGRID_ThlThv, .TRUE. )
+ CALL LES_MEAN_SUBGRID( -XG/PTHVREF/3.*ZA*ZFLX, X_LES_SUBGRID_ThlPz,.TRUE.)
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+!
+!* 8.4 <Rnp Rnp>
+!
+ ! Computes the horizontal variance <Rnp Rnp>
+ IF (.NOT. L2D) THEN
+ ZFLX(:,:,:) = XCHV * PLM(:,:,:) * PLEPS(:,:,:) * &
+ ( GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX)**2 + &
+ GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY)**2 )
+ ELSE
+ ZFLX(:,:,:) = XCHV * PLM(:,:,:) * PLEPS(:,:,:) * &
+ ( GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX)**2 )
+ END IF
+!
+! Compute the flux at the first inner U-point with an uncentred vertical
+! gradient
+ ZFLX(:,:,IKB:IKB) = XCHV * PLM(:,:,IKB:IKB) &
+ * PLEPS(:,:,IKB:IKB) * ( &
+ ( MXF(DXM(PRM(:,:,IKB:IKB,1)) * PINV_PDXX(:,:,IKB:IKB)) &
+ - ( ZCOEFF(:,:,IKB+2:IKB+2)*PRM(:,:,IKB+2:IKB+2,1) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PRM(:,:,IKB+1:IKB+1,1) &
+ +ZCOEFF(:,:,IKB :IKB )*PRM(:,:,IKB :IKB ,1) &
+ ) * 0.5 * ( PDZX(:,:,IKB+1:IKB+1)+PDZX(:,:,IKB:IKB) ) &
+ / MXF(PDXX(:,:,IKB:IKB)) &
+ ) ** 2 + &
+ ( MYF(DYM(PRM(:,:,IKB:IKB,1)) * PINV_PDYY(:,:,IKB:IKB)) &
+ - ( ZCOEFF(:,:,IKB+2:IKB+2)*PRM(:,:,IKB+2:IKB+2,1) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PRM(:,:,IKB+1:IKB+1,1) &
+ +ZCOEFF(:,:,IKB :IKB )*PRM(:,:,IKB :IKB ,1) &
+ ) * 0.5 * ( PDZY(:,:,IKB+1:IKB+1)+PDZY(:,:,IKB:IKB) ) &
+ / MYF(PDYY(:,:,IKB:IKB)) &
+ ) ** 2 )
+!
+ ZFLX(:,:,IKB-1) = ZFLX(:,:,IKB)
+ !
+ IF ( KRRL > 0 ) THEN
+ ZWORK(:,:,:) = ZWORK(:,:,:)+ PAMOIST(:,:,:) * PAMOIST(:,:,:) * ZFLX(:,:,:)
+ END IF
+ !
+ ! stores <Rnp Rnp>
+ IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
+ TZFIELD%CMNHNAME = 'R_HVAR'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'R_HVAR'
+ TZFIELD%CUNITS = 'kg2 kg-2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_R_HVAR'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+ END IF
+ !
+ ! Storage in the LES configuration (addition to TURB_VER computation)
+ !
+ IF (LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( ZFLX, X_LES_SUBGRID_Rt2, .TRUE. )
+ CALL LES_MEAN_SUBGRID( MZF(PWM)*ZFLX, X_LES_RES_W_SBG_Rt2, .TRUE. )
+ CALL LES_MEAN_SUBGRID( ZA*ZFLX, X_LES_SUBGRID_RtThv, .TRUE. )
+ CALL LES_MEAN_SUBGRID( -XG/PTHVREF/3.*ZA*ZFLX, X_LES_SUBGRID_RtPz,.TRUE.)
+ CALL LES_MEAN_SUBGRID( -2.*XCTD*SQRT(PTKEM)*ZFLX/PLEPS, X_LES_SUBGRID_DISS_Rt2, .TRUE. )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+ !
+ END IF
+!
+! 8.5 Complete the Sigma_s computation:
+!
+ IF ( KRRL > 0 ) THEN
+ !
+ PSIGS(:,:,:)=PSIGS(:,:,:)*PSIGS(:,:,:) + ZWORK(:,:,:)
+ ! Extrapolate PSIGS at the ground and at the top
+ PSIGS(:,:,IKB-1) = PSIGS(:,:,IKB)
+ PSIGS(:,:,IKE+1) = PSIGS(:,:,IKE)
+ PSIGS(:,:,:) = SQRT(MAX ( PSIGS(:,:,:),1.E-12) )
+ END IF
+!
+END IF
+!
+!
+!
+END SUBROUTINE TURB_HOR_THERMO_CORR
+END MODULE MODE_TURB_HOR_THERMO_CORR
diff --git a/src/common/turb/mode_turb_hor_thermo_flux.F90 b/src/common/turb/mode_turb_hor_thermo_flux.F90
new file mode 100644
index 0000000000000000000000000000000000000000..97a74596e8bd6ef4086862a7554d27c61c86acf7
--- /dev/null
+++ b/src/common/turb/mode_turb_hor_thermo_flux.F90
@@ -0,0 +1,688 @@
+!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_HOR_THERMO_FLUX
+IMPLICIT NONE
+CONTAINS
+! ################################################################
+ SUBROUTINE TURB_HOR_THERMO_FLUX(KSPLT, KRR, KRRL, KRRI, &
+ OTURB_FLX,OSUBG_COND, &
+ TPFILE, &
+ PK,PINV_PDXX,PINV_PDYY,PINV_PDZZ,PMZM_PRHODJ, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ PDIRCOSXW,PDIRCOSYW, &
+ PRHODJ, &
+ PSFTHM,PSFRM, &
+ PWM,PTHLM,PRM, &
+ PATHETA,PAMOIST,PSRCM,PFRAC_ICE, &
+ PRTHLS,PRRS )
+! ################################################################
+!
+!
+!!**** *TURB_HOR* -routine to compute the source terms in the meso-NH
+!! model equations due to the non-vertical turbulent fluxes.
+!!
+!! PURPOSE
+!! -------
+!!
+!! see TURB_HOR
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!!
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Aug , 1997 (V. Saravane) spliting of TURB_HOR
+!! Nov 27, 1997 (V. Masson) clearing of the routine
+!! Feb. 18, 1998 (J. Stein) bug for v'RC'
+!! Oct 18, 2000 (V. Masson) LES computations + LFLAT switch
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! Feb 20, 2003 (JP Pinty) Add PFRAC_ICE
+!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
+!! change of YCOMMENT
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+USE MODD_CONF
+USE MODD_CTURB
+USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_PARAMETERS
+USE MODD_LES
+!
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
+!
+USE MODI_GRADIENT_M
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_SHUMAN
+USE MODI_LES_MEAN_SUBGRID
+!!USE MODE_EMOIST, ONLY: EMOIST
+!!USE MODE_ETHETA, ONLY: ETHETA
+!
+USE MODI_SECOND_MNH
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declaration of arguments
+!
+!
+!
+INTEGER, INTENT(IN) :: KSPLT ! split process index
+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
+LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for sub-grid
+! condensation
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PK ! Turbulent diffusion doef.
+ ! PK = PLM * SQRT(PTKEM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDXX ! 1./PDXX
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDYY ! 1./PDYY
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDZZ ! 1./PDZZ
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PMZM_PRHODJ ! MZM(PRHODJ)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSXW, PDIRCOSYW
+! Director Cosinus along x, y and z directions at surface w-point
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density * grid volume
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHM,PSFRM
+!
+! Variables at t-1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! mixing ratios at t-1,
+ ! where PRM(:,:,:,1) = conservative mixing ratio
+!
+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(INOUT) :: PRTHLS
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS ! var. at t+1 -split-
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) &
+ :: ZFLX,ZFLXC
+ ! work arrays
+!
+!! REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) :: ZVPTV
+INTEGER :: IKB,IKE,IKU
+ ! Index values for the Beginning and End
+ ! mass points of the domain
+REAL, DIMENSION(SIZE(PDZZ,1),SIZE(PDZZ,2),1+JPVEXT:3+JPVEXT) :: ZCOEFF
+ ! coefficients for the uncentred gradient
+ ! computation near the ground
+!
+REAL :: ZTIME1, ZTIME2
+TYPE(TFIELDDATA) :: TZFIELD
+! ---------------------------------------------------------------------------
+!
+!* 1. PRELIMINARY COMPUTATIONS
+! ------------------------
+!
+IKB = 1+JPVEXT
+IKE = SIZE(PTHLM,3)-JPVEXT
+IKU = SIZE(PTHLM,3)
+!
+!
+! compute the coefficients for the uncentred gradient computation near the
+! ground
+ZCOEFF(:,:,IKB+2)= - PDZZ(:,:,IKB+1) / &
+ ( (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) * PDZZ(:,:,IKB+2) )
+ZCOEFF(:,:,IKB+1)= (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) / &
+ ( PDZZ(:,:,IKB+1) * PDZZ(:,:,IKB+2) )
+ZCOEFF(:,:,IKB)= - (PDZZ(:,:,IKB+2)+2.*PDZZ(:,:,IKB+1)) / &
+ ( (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) * PDZZ(:,:,IKB+1) )
+!
+!* 2. < U' THETA'l >
+! --------------
+!
+!
+ZFLX(:,:,:) = -XCSHF * MXM( PK ) * GX_M_U(1,IKU,1,PTHLM,PDXX,PDZZ,PDZX)
+ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
+!
+! Compute the flux at the first inner U-point with an uncentred vertical
+! gradient
+ZFLX(:,:,IKB:IKB) = -XCSHF * MXM( PK(:,:,IKB:IKB) ) * &
+ ( DXM(PTHLM(:,:,IKB:IKB)) * PINV_PDXX(:,:,IKB:IKB) &
+ -MXM( ZCOEFF(:,:,IKB+2:IKB+2)*PTHLM(:,:,IKB+2:IKB+2) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PTHLM(:,:,IKB+1:IKB+1) &
+ +ZCOEFF(:,:,IKB :IKB )*PTHLM(:,:,IKB :IKB )) &
+ *0.5* ( PDZX(:,:,IKB+1:IKB+1)+PDZX(:,:,IKB:IKB)) &
+ * PINV_PDXX(:,:,IKB:IKB) )
+! extrapolates the flux under the ground so that the vertical average with
+! the IKB flux gives the ground value ( warning the tangential surface
+! flux has been set to 0 for the moment !! to be improved )
+ZFLX(:,:,IKB-1:IKB-1) = 2. * MXM( SPREAD( PSFTHM(:,:)* PDIRCOSXW(:,:), 3,1) ) &
+ - ZFLX(:,:,IKB:IKB)
+!
+! Add this source to the Theta_l sources
+!
+IF (.NOT. LFLAT) THEN
+ PRTHLS(:,:,:) = PRTHLS &
+ - DXF( MXM(PRHODJ) * ZFLX * PINV_PDXX ) &
+ + DZF( PMZM_PRHODJ *MXF(PDZX*(MZM(ZFLX * PINV_PDXX))) * PINV_PDZZ )
+ELSE
+ PRTHLS(:,:,:) = PRTHLS - DXF( MXM(PRHODJ) * ZFLX * PINV_PDXX )
+END IF
+!
+! Compute the equivalent tendancy for Rc and Ri
+!
+IF ( KRRL >= 1 ) THEN
+ IF (.NOT. LFLAT) THEN
+ ZFLXC = 2.*( MXF( MXM( PRHODJ*PATHETA*PSRCM )*ZFLX ) &
+ +MZF( MZM( PRHODJ*PATHETA*PSRCM )*MXF( &
+ PDZX*(MZM( ZFLX*PINV_PDXX )) ) )&
+ )
+ IF ( KRRI >= 1 ) THEN
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
+ (- DXF( MXM( PRHODJ*PATHETA*PSRCM )*ZFLX*PINV_PDXX ) &
+ + DZF( MZM( PRHODJ*PATHETA*PSRCM )*MXF( PDZX*(MZM( ZFLX*PINV_PDXX )) )&
+ *PINV_PDZZ ) &
+ )*(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,4) = PRRS(:,:,:,4) + 2. * &
+ (- DXF( MXM( PRHODJ*PATHETA*PSRCM )*ZFLX*PINV_PDXX ) &
+ + DZF( MZM( PRHODJ*PATHETA*PSRCM )*MXF( PDZX*(MZM( ZFLX*PINV_PDXX )) )&
+ *PINV_PDZZ ) &
+ )*PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
+ (- DXF( MXM( PRHODJ*PATHETA*PSRCM )*ZFLX*PINV_PDXX ) &
+ + DZF( MZM( PRHODJ*PATHETA*PSRCM )*MXF( PDZX*(MZM( ZFLX*PINV_PDXX )) )&
+ *PINV_PDZZ ) &
+ )
+ END IF
+ ELSE
+ ZFLXC = 2.*MXF( MXM( PRHODJ*PATHETA*PSRCM )*ZFLX )
+ IF ( KRRI >= 1 ) THEN
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
+ DXF( MXM( PRHODJ*PATHETA*PSRCM )*ZFLX*PINV_PDXX )*(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * &
+ DXF( MXM( PRHODJ*PATHETA*PSRCM )*ZFLX*PINV_PDXX )*PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
+ DXF( MXM( PRHODJ*PATHETA*PSRCM )*ZFLX*PINV_PDXX )
+ END IF
+ END IF
+END IF
+!
+!! stores this flux in ZWORK to compute later <U' VPT'>
+!!ZWORK(:,:,:) = ZFLX(:,:,:)
+!
+! stores the horizontal <U THl>
+IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+ TZFIELD%CMNHNAME = 'UTHL_FLX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'UTHL_FLX'
+ TZFIELD%CUNITS = 'K m s-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_UTHL_FLX'
+ TZFIELD%NGRID = 2
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+END IF
+!
+IF (KSPLT==1 .AND. LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( MXF(ZFLX), X_LES_SUBGRID_UThl )
+ CALL LES_MEAN_SUBGRID( MZF(MXF(GX_W_UW(PWM,PDXX,PDZZ,PDZX)*MZM(ZFLX))),&
+ X_LES_RES_ddxa_W_SBG_UaThl , .TRUE. )
+ CALL LES_MEAN_SUBGRID( GX_M_M(PTHLM,PDXX,PDZZ,PDZX)*MXF(ZFLX),&
+ X_LES_RES_ddxa_Thl_SBG_UaThl , .TRUE. )
+ IF (KRR>=1) THEN
+ CALL LES_MEAN_SUBGRID( GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX)*MXF(ZFLX), &
+ X_LES_RES_ddxa_Rt_SBG_UaThl , .TRUE. )
+ END IF
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+!
+!* 3. < U' R'np >
+! -----------
+IF (KRR/=0) THEN
+ !
+ ZFLX(:,:,:) = -XCHF * MXM( PK ) * GX_M_U(1,IKU,1,PRM(:,:,:,1),PDXX,PDZZ,PDZX)
+ ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
+!
+! Compute the flux at the first inner U-point with an uncentred vertical
+! gradient
+ ZFLX(:,:,IKB:IKB) = -XCHF * MXM( PK(:,:,IKB:IKB) ) * &
+ ( DXM(PRM(:,:,IKB:IKB,1)) * PINV_PDXX(:,:,IKB:IKB) &
+ -MXM( ZCOEFF(:,:,IKB+2:IKB+2)*PRM(:,:,IKB+2:IKB+2,1) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PRM(:,:,IKB+1:IKB+1,1) &
+ +ZCOEFF(:,:,IKB :IKB )*PRM(:,:,IKB :IKB ,1)) &
+ *0.5* ( PDZX(:,:,IKB+1:IKB+1)+PDZX(:,:,IKB:IKB)) &
+ * PINV_PDXX(:,:,IKB:IKB) )
+! extrapolates the flux under the ground so that the vertical average with
+! the IKB flux gives the ground value ( warning the tangential surface
+! flux has been set to 0 for the moment !! to be improved )
+ ZFLX(:,:,IKB-1:IKB-1) = 2. * MXM( SPREAD( PSFRM(:,:)* PDIRCOSXW(:,:), 3,1) ) &
+ - ZFLX(:,:,IKB:IKB)
+ !
+ ! Add this source to the conservative mixing ratio sources
+ !
+ IF (.NOT. LFLAT) THEN
+ PRRS(:,:,:,1) = PRRS(:,:,:,1) &
+ - DXF( MXM(PRHODJ) * ZFLX * PINV_PDXX ) &
+ + DZF( PMZM_PRHODJ *MXF(PDZX*(MZM(ZFLX * PINV_PDXX))) * PINV_PDZZ )
+ ELSE
+ PRRS(:,:,:,1) = PRRS(:,:,:,1) - DXF( MXM(PRHODJ) * ZFLX * PINV_PDXX )
+ END IF
+ !
+ ! Compute the equivalent tendancy for Rc and Ri
+ !
+ IF ( KRRL >= 1 ) THEN
+ IF (.NOT. LFLAT) THEN
+ ZFLXC = ZFLXC &
+ + 2.*( MXF( MXM( PRHODJ*PAMOIST*PSRCM )*ZFLX ) &
+ +MZF( MZM( PRHODJ*PAMOIST*PSRCM )*MXF( &
+ PDZX*(MZM( ZFLX*PINV_PDXX )) ) )&
+ )
+ IF ( KRRI >= 1 ) THEN
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
+ (- DXF( MXM( PRHODJ*PAMOIST*PSRCM )*ZFLX*PINV_PDXX ) &
+ + DZF( MZM( PRHODJ*PAMOIST*PSRCM )*MXF( PDZX*(MZM( ZFLX*PINV_PDXX )) )&
+ *PINV_PDZZ ) &
+ )*(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
+ (- DXF( MXM( PRHODJ*PAMOIST*PSRCM )*ZFLX*PINV_PDXX ) &
+ + DZF( MZM( PRHODJ*PAMOIST*PSRCM )*MXF( PDZX*(MZM( ZFLX*PINV_PDXX )) )&
+ *PINV_PDZZ ) &
+ )*PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
+ (- DXF( MXM( PRHODJ*PAMOIST*PSRCM )*ZFLX*PINV_PDXX ) &
+ + DZF( MZM( PRHODJ*PAMOIST*PSRCM )*MXF( PDZX*(MZM( ZFLX*PINV_PDXX )) )&
+ *PINV_PDZZ ) &
+ )
+ END IF
+ ELSE
+ ZFLXC = ZFLXC + 2.*MXF( MXM( PRHODJ*PAMOIST*PSRCM )*ZFLX )
+ IF ( KRRI >= 1 ) THEN
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
+ DXF( MXM( PRHODJ*PAMOIST*PSRCM )*ZFLX*PINV_PDXX )*(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * &
+ DXF( MXM( PRHODJ*PAMOIST*PSRCM )*ZFLX*PINV_PDXX )*PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
+ DXF( MXM( PRHODJ*PAMOIST*PSRCM )*ZFLX*PINV_PDXX )
+ END IF
+ END IF
+ END IF
+ !
+ ! stores the horizontal <U Rnp>
+ IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+ TZFIELD%CMNHNAME = 'UR_FLX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'UR_FLX'
+ TZFIELD%CUNITS = 'kg kg-1 m s-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_UR_FLX'
+ TZFIELD%NGRID = 2
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+ END IF
+ !
+ IF (KSPLT==1 .AND. LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( MXF(ZFLX), X_LES_SUBGRID_URt )
+ CALL LES_MEAN_SUBGRID( MZF(MXF(GX_W_UW(PWM,PDXX,PDZZ,PDZX)*MZM(ZFLX))),&
+ X_LES_RES_ddxa_W_SBG_UaRt , .TRUE. )
+ CALL LES_MEAN_SUBGRID( GX_M_M(PTHLM,PDXX,PDZZ,PDZX)*MXF(ZFLX),&
+ X_LES_RES_ddxa_Thl_SBG_UaRt , .TRUE. )
+ CALL LES_MEAN_SUBGRID( GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX)*MXF(ZFLX),&
+ X_LES_RES_ddxa_Rt_SBG_UaRt , .TRUE. )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+!
+ !
+ IF (KRRL>0 .AND. KSPLT==1 .AND. LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID(MXF(ZFLXC), X_LES_SUBGRID_URc )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+!
+END IF
+!
+!* 4. < U' TPV' >
+! -----------
+!
+!! to be tested later
+!!IF (KRR/=0) THEN
+!! ! here ZFLX= <U'Rnp'> and ZWORK= <U'Thetal'>
+!! !
+!! ZVPTU(:,:,:) = &
+!! ZWORK(:,:,:)*MXM(ETHETA(KRR,KRRI,PTHLT,PEXNREF,PRT,PLOCPT,PSRCM)) + &
+!! ZFLX(:,:,:)*MXM(EMOIST(KRR,KRRI,PTHLT,PEXNREF,PRT,PLOCPT,PSRCM))
+!! !
+!! ! stores the horizontal <U VPT>
+!! IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+!! TZFIELD%CMNHNAME = 'UVPT_FLX'
+!! TZFIELD%CSTDNAME = ''
+!! TZFIELD%CLONGNAME = 'UVPT_FLX'
+!! TZFIELD%CUNITS = 'K m s-1'
+!! TZFIELD%CDIR = 'XY'
+!! TZFIELD%CCOMMENT = 'X_Y_Z_UVPT_FLX'
+!! TZFIELD%NGRID = 2
+!! TZFIELD%NTYPE = TYPEREAL
+!! TZFIELD%NDIMS = 3
+!! TZFIELD%LTIMEDEP = .TRUE.
+!! CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZVPTU)
+!! END IF
+!!!
+!!ELSE
+!! ZVPTU(:,:,:)=ZWORK(:,:,:)
+!!END IF
+!
+!
+!* 5. < V' THETA'l >
+! --------------
+!
+!
+IF (.NOT. L2D) THEN
+ ZFLX(:,:,:) = -XCSHF * MYM( PK ) * GY_M_V(1,IKU,1,PTHLM,PDYY,PDZZ,PDZY)
+ ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
+ELSE
+ ZFLX(:,:,:) = 0.
+END IF
+!
+!
+! Compute the flux at the first inner U-point with an uncentred vertical
+! gradient
+ZFLX(:,:,IKB:IKB) = -XCSHF * MYM( PK(:,:,IKB:IKB) ) * &
+ ( DYM(PTHLM(:,:,IKB:IKB)) * PINV_PDYY(:,:,IKB:IKB) &
+ -MYM( ZCOEFF(:,:,IKB+2:IKB+2)*PTHLM(:,:,IKB+2:IKB+2) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PTHLM(:,:,IKB+1:IKB+1) &
+ +ZCOEFF(:,:,IKB :IKB )*PTHLM(:,:,IKB :IKB ) ) &
+ *0.5* ( PDZY(:,:,IKB+1:IKB+1)+PDZY(:,:,IKB:IKB)) &
+ * PINV_PDYY(:,:,IKB:IKB) )
+! extrapolates the flux under the ground so that the vertical average with
+! the IKB flux gives the ground value ( warning the tangential surface
+! flux has been set to 0 for the moment !! to be improved )
+ZFLX(:,:,IKB-1:IKB-1) = 2. * MYM( SPREAD( PSFTHM(:,:)* PDIRCOSYW(:,:), 3,1) ) &
+ - ZFLX(:,:,IKB:IKB)
+!
+! Add this source to the Theta_l sources
+!
+IF (.NOT. L2D) THEN
+ IF (.NOT. LFLAT) THEN
+ PRTHLS(:,:,:) = PRTHLS &
+ - DYF( MYM(PRHODJ) * ZFLX * PINV_PDYY ) &
+ + DZF( PMZM_PRHODJ *MYF(PDZY*(MZM(ZFLX * PINV_PDYY))) * PINV_PDZZ )
+ ELSE
+ PRTHLS(:,:,:) = PRTHLS - DYF( MYM(PRHODJ) * ZFLX * PINV_PDYY )
+ END IF
+END IF
+!
+! Compute the equivalent tendancy for Rc and Ri
+!
+!IF ( OSUBG_COND .AND. KRRL > 0 .AND. .NOT. L2D) THEN
+IF ( KRRL >= 1 .AND. .NOT. L2D) THEN
+ IF (.NOT. LFLAT) THEN
+ ZFLXC = 2.*( MYF( MYM( PRHODJ*PATHETA*PSRCM )*ZFLX ) &
+ +MZF( MZM( PRHODJ*PATHETA*PSRCM )*MYF( &
+ PDZY*(MZM( ZFLX*PINV_PDYY )) ) )&
+ )
+ IF ( KRRI >= 1 ) THEN
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
+ (- DYF( MYM( PRHODJ*PATHETA*PSRCM )*ZFLX*PINV_PDYY ) &
+ + DZF( MZM( PRHODJ*PATHETA*PSRCM )*MYF( PDZY*(MZM( ZFLX*PINV_PDYY )) )&
+ *PINV_PDZZ ) &
+ )*(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,4) = PRRS(:,:,:,4) + 2. * &
+ (- DYF( MYM( PRHODJ*PATHETA*PSRCM )*ZFLX*PINV_PDYY ) &
+ + DZF( MZM( PRHODJ*PATHETA*PSRCM )*MYF( PDZY*(MZM( ZFLX*PINV_PDYY )) )&
+ *PINV_PDZZ ) &
+ )*PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
+ (- DYF( MYM( PRHODJ*PATHETA*PSRCM )*ZFLX*PINV_PDYY ) &
+ + DZF( MZM( PRHODJ*PATHETA*PSRCM )*MYF( PDZY*(MZM( ZFLX*PINV_PDYY )) )&
+ *PINV_PDZZ ) &
+ )
+ END IF
+ ELSE
+ ZFLXC = 2.*MYF( MYM( PRHODJ*PATHETA*PSRCM )*ZFLX )
+ IF ( KRRI >= 1 ) THEN
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
+ DYF( MYM( PRHODJ*PATHETA*PSRCM )*ZFLX*PINV_PDYY )*(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * &
+ DYF( MYM( PRHODJ*PATHETA*PSRCM )*ZFLX*PINV_PDYY )*PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
+ DYF( MYM( PRHODJ*PATHETA*PSRCM )*ZFLX*PINV_PDYY )
+ END IF
+ END IF
+END IF
+!
+!! stores this flux in ZWORK to compute later <V' VPT'>
+!!ZWORK(:,:,:) = ZFLX(:,:,:)
+!
+! stores the horizontal <V THl>
+IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+ TZFIELD%CMNHNAME = 'VTHL_FLX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'VTHL_FLX'
+ TZFIELD%CUNITS = 'K m s-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_VTHL_FLX'
+ TZFIELD%NGRID = 3
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+END IF
+!
+IF (KSPLT==1 .AND. LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( MYF(ZFLX), X_LES_SUBGRID_VThl )
+ CALL LES_MEAN_SUBGRID( MZF(MYF(GY_W_VW(PWM,PDYY,PDZZ,PDZY)*MZM(ZFLX))),&
+ X_LES_RES_ddxa_W_SBG_UaThl , .TRUE. )
+ CALL LES_MEAN_SUBGRID( GY_M_M(PTHLM,PDYY,PDZZ,PDZY)*MYF(ZFLX),&
+ X_LES_RES_ddxa_Thl_SBG_UaThl , .TRUE. )
+ IF (KRR>=1) THEN
+ CALL LES_MEAN_SUBGRID( GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY)*MYF(ZFLX),&
+ X_LES_RES_ddxa_Rt_SBG_UaThl , .TRUE. )
+ END IF
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+!
+!
+!* 6. < V' R'np >
+! -----------
+!
+IF (KRR/=0) THEN
+ !
+ IF (.NOT. L2D) THEN
+ ZFLX(:,:,:) = -XCHF * MYM( PK ) * GY_M_V(1,IKU,1,PRM(:,:,:,1),PDYY,PDZZ,PDZY)
+ ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
+ ELSE
+ ZFLX(:,:,:) = 0.
+ END IF
+!
+! Compute the flux at the first inner U-point with an uncentred vertical
+! gradient
+ ZFLX(:,:,IKB:IKB) = -XCHF * MYM( PK(:,:,IKB:IKB) ) * &
+ ( DYM(PRM(:,:,IKB:IKB,1)) * PINV_PDYY(:,:,IKB:IKB) &
+ -MYM( ZCOEFF(:,:,IKB+2:IKB+2)*PRM(:,:,IKB+2:IKB+2,1) &
+ +ZCOEFF(:,:,IKB+1:IKB+1)*PRM(:,:,IKB+1:IKB+1,1) &
+ +ZCOEFF(:,:,IKB :IKB )*PRM(:,:,IKB :IKB ,1) ) &
+ *0.5* ( PDZY(:,:,IKB+1:IKB+1)+PDZY(:,:,IKB:IKB)) &
+ * PINV_PDYY(:,:,IKB:IKB) )
+! extrapolates the flux under the ground so that the vertical average with
+! the IKB flux gives the ground value ( warning the tangential surface
+! flux has been set to 0 for the moment !! to be improved )
+ ZFLX(:,:,IKB-1:IKB-1) = 2. * MYM( SPREAD( PSFRM(:,:)* PDIRCOSYW(:,:), 3,1) ) &
+ - ZFLX(:,:,IKB:IKB)
+ !
+ ! Add this source to the conservative mixing ratio sources
+ !
+ IF (.NOT. L2D) THEN
+ IF (.NOT. LFLAT) THEN
+ PRRS(:,:,:,1) = PRRS(:,:,:,1) &
+ - DYF( MYM(PRHODJ) * ZFLX * PINV_PDYY ) &
+
+ + DZF( PMZM_PRHODJ *MYF(PDZY*(MZM(ZFLX * PINV_PDYY))) * PINV_PDZZ )
+ ELSE
+ PRRS(:,:,:,1) = PRRS(:,:,:,1) - DYF( MYM(PRHODJ) * ZFLX * PINV_PDYY )
+ END IF
+ END IF
+ !
+ ! Compute the equivalent tendancy for Rc and Ri
+ !
+ IF ( KRRL >= 1 .AND. .NOT. L2D) THEN ! Sub-grid condensation
+ IF (.NOT. LFLAT) THEN
+ ZFLXC = ZFLXC &
+ + 2.*( MXF( MYM( PRHODJ*PAMOIST*PSRCM )*ZFLX ) &
+ + MZF( MZM( PRHODJ*PAMOIST*PSRCM )*MYF( &
+ PDZY*(MZM( ZFLX*PINV_PDYY )) ) )&
+ )
+ IF ( KRRI >= 1 ) THEN
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
+ (- DYF( MYM( PRHODJ*PAMOIST*PSRCM )*ZFLX/PDYY ) &
+ + DZF( MZM( PRHODJ*PAMOIST*PSRCM )*MYF( PDZY*(MZM( ZFLX*PINV_PDYY )) )&
+ * PINV_PDZZ ) &
+ )*(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,4) = PRRS(:,:,:,4) + 2. * &
+ (- DYF( MYM( PRHODJ*PAMOIST*PSRCM )*ZFLX/PDYY ) &
+ + DZF( MZM( PRHODJ*PAMOIST*PSRCM )*MYF( PDZY*(MZM( ZFLX*PINV_PDYY )) )&
+ * PINV_PDZZ ) &
+ )*PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
+ (- DYF( MYM( PRHODJ*PAMOIST*PSRCM )*ZFLX/PDYY ) &
+ + DZF( MZM( PRHODJ*PAMOIST*PSRCM )*MYF( PDZY*(MZM( ZFLX*PINV_PDYY )) )&
+ * PINV_PDZZ ) &
+ )
+ END IF
+ ELSE
+ ZFLXC = ZFLXC + 2.*MXF( MYM( PRHODJ*PAMOIST*PSRCM )*ZFLX )
+ IF ( KRRI >= 1 ) THEN
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
+ DYF( MYM( PRHODJ*PAMOIST*PSRCM )*ZFLX/PDYY )*(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * &
+ DYF( MYM( PRHODJ*PAMOIST*PSRCM )*ZFLX/PDYY )*PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
+ DYF( MYM( PRHODJ*PAMOIST*PSRCM )*ZFLX/PDYY )
+ END IF
+ END IF
+ END IF
+ !
+ ! stores the horizontal <V Rnp>
+ IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+ TZFIELD%CMNHNAME = 'VR_FLX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'VR_FLX'
+ TZFIELD%CUNITS = 'kg kg-1 m s-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_VR_FLX'
+ TZFIELD%NGRID = 3
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+ END IF
+ !
+ IF (KSPLT==1 .AND. LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( MYF(ZFLX), X_LES_SUBGRID_VRt )
+ CALL LES_MEAN_SUBGRID( MZF(MYF(GY_W_VW(PWM,PDYY,PDZZ,PDZY)*MZM(ZFLX))),&
+ X_LES_RES_ddxa_W_SBG_UaRt , .TRUE. )
+ CALL LES_MEAN_SUBGRID( GY_M_M(PTHLM,PDYY,PDZZ,PDZY)*MYF(ZFLX), &
+ X_LES_RES_ddxa_Thl_SBG_UaRt , .TRUE. )
+ CALL LES_MEAN_SUBGRID( GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY)*MYF(ZFLX), &
+ X_LES_RES_ddxa_Rt_SBG_UaRt , .TRUE. )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+!
+ !
+ IF (KRRL>0 .AND. KSPLT==1 .AND. LLES_CALL) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID(MYF(ZFLXC), X_LES_SUBGRID_VRc )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+ !
+END IF
+!
+!* 7. < V' TPV' >
+! -----------
+!
+!! to be tested later
+!!IF (KRR/=0) THEN
+!! ! here ZFLX= <V'R'np> and ZWORK= <V'Theta'l>
+!! !
+!! IF (.NOT. L2D) THEN &
+!! ZVPTV(:,:,:) = &
+!! ZWORK(:,:,:)*MYM(ETHETA(KRR,KRRI,PTHLT,PEXNREF,PRT,PLOCPT,PSRCM)) + &
+!! ZFLX(:,:,:)*MYM(EMOIST(KRR,KRRI,PTHLT,PEXNREF,PRT,PLOCPT,PSRCM))
+!! ELSE
+!! ZVPTV(:,:,:) = 0.
+!! END IF
+!! !
+!! ! stores the horizontal <V VPT>
+!! IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+!! TZFIELD%CMNHNAME = 'VVPT_FLX'
+!! TZFIELD%CSTDNAME = ''
+!! TZFIELD%CLONGNAME = 'VVPT_FLX'
+!! TZFIELD%CUNITS = 'K m s-1'
+!! TZFIELD%CDIR = 'XY'
+!! TZFIELD%CCOMMENT = 'X_Y_Z_VVPT_FLX'
+!! TZFIELD%NGRID = 3
+!! TZFIELD%NTYPE = TYPEREAL
+!! TZFIELD%NDIMS = 3
+!! TZFIELD%LTIMEDEP = .TRUE.
+!! CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZVPTV)
+!! END IF
+!!!
+!!ELSE
+!! ZVPTV(:,:,:)=ZWORK(:,:,:)
+!!END IF
+!
+!
+END SUBROUTINE TURB_HOR_THERMO_FLUX
+END MODULE MODE_TURB_HOR_THERMO_FLUX
diff --git a/src/common/turb/mode_turb_hor_tke.F90 b/src/common/turb/mode_turb_hor_tke.F90
new file mode 100644
index 0000000000000000000000000000000000000000..5ff7a0029077c979ca71d3ed008b7d165cc340f8
--- /dev/null
+++ b/src/common/turb/mode_turb_hor_tke.F90
@@ -0,0 +1,213 @@
+!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 MODE_TURB_HOR_TKE
+IMPLICIT NONE
+CONTAINS
+ SUBROUTINE TURB_HOR_TKE(KSPLT, &
+ PDXX, PDYY, PDZZ,PDZX,PDZY, &
+ PINV_PDXX, PINV_PDYY, PINV_PDZZ, PMZM_PRHODJ, &
+ PK, PRHODJ, PTKEM, &
+ PTRH )
+! ################################################################
+!
+!
+!!**** *TURB_HOR_TKE* computes the horizontal turbulant transports of Tke
+!!
+!! PURPOSE
+!! -------
+
+!!** METHOD
+!! ------
+!!
+!!
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original Aug 29, 1994
+!! Mar 07 2001 (V. Masson and J. Stein) new routine
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CONF
+USE MODD_CST
+USE MODD_CTURB
+USE MODD_PARAMETERS
+USE MODD_LES
+!
+!
+USE MODI_SHUMAN
+USE MODI_GRADIENT_M
+USE MODI_LES_MEAN_SUBGRID
+!
+USE MODI_SECOND_MNH
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declaration of arguments
+!
+!
+INTEGER, INTENT(IN) :: KSPLT ! current split index
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PK ! Turbulent diffusion doef.
+ ! PK = PLM * SQRT(PTKEM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDXX ! 1./PDXX
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDYY ! 1./PDYY
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDZZ ! 1./PDZZ
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PMZM_PRHODJ ! MZM(PRHODJ)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density * grid volume
+!
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at time t- dt
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PTRH ! horizontal transport of Tke
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+INTEGER :: IKB, IKU
+!
+REAL, DIMENSION(SIZE(PDZZ,1),SIZE(PDZZ,2),1+JPVEXT:3+JPVEXT) :: ZCOEFF
+ ! coefficients for the uncentred gradient
+ ! computation near the ground
+!
+REAL, DIMENSION(SIZE(PTKEM,1),SIZE(PTKEM,2),SIZE(PTKEM,3)):: ZFLX
+!
+REAL :: ZTIME1, ZTIME2
+! ---------------------------------------------------------------------------
+!
+!* 1. PRELIMINARY COMPUTATIONS
+! ------------------------
+!
+IKB = 1.+JPVEXT
+IKU = SIZE(PTKEM,3)
+!
+! compute the coefficients for the uncentred gradient computation near the
+! ground
+!
+ZCOEFF(:,:,IKB+2)= - PDZZ(:,:,IKB+1) / &
+ ( (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) * PDZZ(:,:,IKB+2) )
+ZCOEFF(:,:,IKB+1)= (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) / &
+ ( PDZZ(:,:,IKB+1) * PDZZ(:,:,IKB+2) )
+ZCOEFF(:,:,IKB)= - (PDZZ(:,:,IKB+2)+2.*PDZZ(:,:,IKB+1)) / &
+ ( (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) * PDZZ(:,:,IKB+1) )
+!
+!--------------------------------------------------------------------
+!
+!* 2. horizontal transport of Tke u'e
+! -------------------------------
+!
+!
+ZFLX = -XCET * MXM(PK) * GX_M_U(1,IKU,1,PTKEM,PDXX,PDZZ,PDZX) ! < u'e >
+!
+! special case near the ground ( uncentred gradient )
+!
+ZFLX(:,:,IKB) = ZCOEFF(:,:,IKB+2)*PTKEM(:,:,IKB+2) &
+ + ZCOEFF(:,:,IKB+1)*PTKEM(:,:,IKB+1) &
+ + ZCOEFF(:,:,IKB )*PTKEM(:,:,IKB )
+!
+ZFLX(:,:,IKB:IKB) = &
+ - XCET * MXM( PK(:,:,IKB:IKB) ) * ( &
+ DXM ( PTKEM(:,:,IKB:IKB) ) * PINV_PDXX(:,:,IKB:IKB) &
+ -MXM ( ZFLX (:,:,IKB:IKB) ) * PINV_PDXX(:,:,IKB:IKB) &
+ * 0.5 * ( PDZX(:,:,IKB+1:IKB+1) + PDZX(:,:,IKB:IKB) ) )
+!
+! extrapolate the fluxes to obtain < u'e > = 0 at the ground
+!
+ZFLX(:,:,IKB-1) = - ZFLX(:,:,IKB)
+!
+! let the same flux at IKU-1 and IKU level
+!
+ZFLX(:,:,IKU) = ZFLX(:,:,IKU-1)
+!
+IF (.NOT. LFLAT) THEN
+ PTRH =-( DXF( MXM(PRHODJ) * ZFLX * PINV_PDXX)&
+ - DZF( PMZM_PRHODJ * MXF( PDZX * MZM(ZFLX*PINV_PDXX)) * PINV_PDZZ)&
+ ) /PRHODJ
+ELSE
+ PTRH =-( DXF( MXM(PRHODJ) * ZFLX * PINV_PDXX)&
+ ) /PRHODJ
+END IF
+!
+IF (LLES_CALL .AND. KSPLT==1) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( MXF(ZFLX), X_LES_SUBGRID_UTke )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+!
+!
+!--------------------------------------------------------------------
+!
+!* 3. horizontal transport of Tke v'e
+! -------------------------------
+!
+IF (.NOT. L2D) THEN
+ ZFLX =-XCET * MYM(PK) * GY_M_V(1,IKU,1,PTKEM,PDYY,PDZZ,PDZY) ! < v'e >
+!
+! special case near the ground ( uncentred gradient )
+!
+ ZFLX(:,:,IKB) = ZCOEFF(:,:,IKB+2)*PTKEM(:,:,IKB+2) &
+ + ZCOEFF(:,:,IKB+1)*PTKEM(:,:,IKB+1) &
+ + ZCOEFF(:,:,IKB )*PTKEM(:,:,IKB )
+!
+ ZFLX(:,:,IKB:IKB) = &
+ - XCET * MYM( PK(:,:,IKB:IKB) ) * ( &
+ DYM ( PTKEM(:,:,IKB:IKB) ) * PINV_PDYY(:,:,IKB:IKB) &
+ - MYM ( ZFLX (:,:,IKB:IKB) ) * PINV_PDYY(:,:,IKB:IKB) &
+ * 0.5 * ( PDZY(:,:,IKB+1:IKB+1) + PDZY(:,:,IKB:IKB) ) )
+!
+! extrapolate the fluxes to obtain < v'e > = 0 at the ground
+!
+ ZFLX(:,:,IKB-1) = - ZFLX(:,:,IKB)
+!
+! let the same flux at IKU-1 and IKU level
+!
+ ZFLX(:,:,IKU) = ZFLX(:,:,IKU-1)
+!
+! complete the explicit turbulent transport
+!
+ IF (.NOT. LFLAT) THEN
+ PTRH = PTRH - ( DYF( MYM(PRHODJ) * ZFLX * PINV_PDYY ) &
+ - DZF( PMZM_PRHODJ * MYF( PDZY * MZM(ZFLX*PINV_PDYY) ) * PINV_PDZZ ) &
+ ) /PRHODJ
+ ELSE
+ PTRH = PTRH - ( DYF( MYM(PRHODJ) * ZFLX * PINV_PDYY ) &
+ ) /PRHODJ
+ END IF
+!
+ IF (LLES_CALL .AND. KSPLT==1) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( MYF(ZFLX), X_LES_SUBGRID_VTke )
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+ END IF
+!
+END IF
+!
+!----------------------------------------------------------------------------
+!
+END SUBROUTINE TURB_HOR_TKE
+END MODULE MODE_TURB_HOR_TKE
diff --git a/src/common/turb/mode_turb_hor_uv.F90 b/src/common/turb/mode_turb_hor_uv.F90
new file mode 100644
index 0000000000000000000000000000000000000000..611679c73c5fa03ddca83a2daddd5cd1707a4d28
--- /dev/null
+++ b/src/common/turb/mode_turb_hor_uv.F90
@@ -0,0 +1,292 @@
+!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_HOR_UV
+IMPLICIT NONE
+CONTAINS
+! ################################################################
+ SUBROUTINE TURB_HOR_UV(KSPLT, &
+ OTURB_FLX, &
+ TPFILE, &
+ PK,PINV_PDXX,PINV_PDYY,PINV_PDZZ,PMZM_PRHODJ, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY, &
+ PDIRCOSZW, &
+ PCOSSLOPE,PSINSLOPE, &
+ PRHODJ, &
+ PCDUEFF,PTAU11M,PTAU12M,PTAU22M,PTAU33M, &
+ PUM,PVM,PUSLOPEM,PVSLOPEM, &
+ PDP, &
+ PRUS,PRVS )
+! ################################################################
+!
+!
+!!**** *TURB_HOR* -routine to compute the source terms in the meso-NH
+!! model equations due to the non-vertical turbulent fluxes.
+!!
+!! PURPOSE
+!! -------
+!!
+!! see TURB_HOR
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!!
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Aug , 1997 (V. Saravane) spliting of TURB_HOR
+!! Nov 27, 1997 (V. Masson) clearing of the routine
+!! Oct 18, 2000 (V. Masson) LES computations + LFLAT switch
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+USE MODD_CONF
+USE MODD_CTURB
+USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_PARAMETERS
+USE MODD_LES
+!
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
+!
+USE MODI_GRADIENT_M
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_SHUMAN
+USE MODE_COEFJ, ONLY: COEFJ
+USE MODI_LES_MEAN_SUBGRID
+!
+USE MODI_SECOND_MNH
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declaration of arguments
+!
+!
+!
+INTEGER, INTENT(IN) :: KSPLT ! split process index
+LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
+ ! turbulent fluxes in the syncronous FM-file
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PK ! Turbulent diffusion doef.
+ ! PK = PLM * SQRT(PTKEM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDXX ! 1./PDXX
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDYY ! 1./PDYY
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDZZ ! 1./PDZZ
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PMZM_PRHODJ ! MZM(PRHODJ)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSZW
+! Director Cosinus along z directions at surface w-point
+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 ! density * grid volume
+!
+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) :: PTAU22M ! <vv> in the same axes
+REAL, DIMENSION(:,:), INTENT(IN) :: PTAU33M ! <ww> in the same axes
+!
+! Variables at t-1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUM,PVM
+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(INOUT) :: PRUS, PRVS ! var. at t+1 -split-
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDP ! TKE production terms
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),SIZE(PUM,3)) &
+ :: ZFLX,ZWORK
+ ! work arrays
+!
+REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2)) ::ZDIRSINZW
+ ! sinus of the angle between the vertical and the normal to the orography
+INTEGER :: IKB,IKE
+ ! Index values for the Beginning and End
+ ! mass points of the domain
+!
+REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),SIZE(PUM,3)) :: GY_U_UV_PUM
+REAL, DIMENSION(SIZE(PVM,1),SIZE(PVM,2),SIZE(PVM,3)) :: GX_V_UV_PVM
+!
+REAL :: ZTIME1, ZTIME2
+TYPE(TFIELDDATA) :: TZFIELD
+! ---------------------------------------------------------------------------
+!
+!* 1. PRELIMINARY COMPUTATIONS
+! ------------------------
+!
+IKB = 1+JPVEXT
+IKE = SIZE(PUM,3)-JPVEXT
+!
+ZDIRSINZW(:,:) = SQRT( 1. - PDIRCOSZW(:,:)**2 )
+!
+GX_V_UV_PVM = GX_V_UV(PVM,PDXX,PDZZ,PDZX)
+IF (.NOT. L2D) GY_U_UV_PUM = GY_U_UV(PUM,PDYY,PDZZ,PDZY)
+!
+!
+!* 12. < U'V'>
+! -------
+!
+!
+IF (.NOT. L2D) THEN
+ ZFLX(:,:,:)= - XCMFS * MYM(MXM(PK)) * &
+ (GY_U_UV_PUM + GX_V_UV_PVM)
+ELSE
+ ZFLX(:,:,:)= - XCMFS * MYM(MXM(PK)) * &
+ (GX_V_UV_PVM)
+END IF
+!
+ZFLX(:,:,IKE+1)= ZFLX(:,:,IKE)
+!
+!
+! Compute the correlation at the first physical level with the following
+! hypothesis du/dz vary in 1/z and w=0 at the ground
+ZFLX(:,:,IKB:IKB) = - XCMFS * MYM(MXM(PK(:,:,IKB:IKB))) * ( &
+ ( DYM( PUM(:,:,IKB:IKB) ) &
+ -MYM( (PUM(:,:,IKB+1:IKB+1)-PUM(:,:,IKB:IKB)) &
+ *(1./MXM(PDZZ(:,:,IKB+1:IKB+1))+1./MXM(PDZZ(:,:,IKB:IKB))))&
+ *0.5*MXM((PDZY(:,:,IKB+1:IKB+1)+PDZY(:,:,IKB:IKB))) &
+ ) / MXM(PDYY(:,:,IKB:IKB)) &
+ +( DXM( PVM(:,:,IKB:IKB) ) &
+ -MXM( (PVM(:,:,IKB+1:IKB+1)-PVM(:,:,IKB:IKB)) &
+ *(1./MYM(PDZZ(:,:,IKB+1:IKB+1))+1./MYM(PDZZ(:,:,IKB:IKB))))&
+ *0.5*MYM((PDZX(:,:,IKB+1:IKB+1)+PDZX(:,:,IKB:IKB))) &
+ ) / MYM(PDXX(:,:,IKB:IKB)) )
+!
+! extrapolates this flux under the ground with the surface flux
+ZFLX(:,:,IKB-1) = &
+ PTAU11M(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * PDIRCOSZW(:,:)**2 &
+ +PTAU12M(:,:) * (PCOSSLOPE(:,:)**2 - PSINSLOPE(:,:)**2) * &
+ PDIRCOSZW(:,:)**2 &
+ -PTAU22M(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) &
+ +PTAU33M(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:)**2 &
+ -PCDUEFF(:,:) * ( &
+ 2. * PUSLOPEM(:,:) * PCOSSLOPE(:,:) * PSINSLOPE(:,:) * &
+ PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
+ +PVSLOPEM(:,:) * (PCOSSLOPE(:,:)**2 - PSINSLOPE(:,:)**2) * ZDIRSINZW(:,:) &
+ )
+!
+ZFLX(:,:,IKB-1:IKB-1) = 2. * MXM( MYM( ZFLX(:,:,IKB-1:IKB-1) ) ) &
+ - ZFLX(:,:,IKB:IKB)
+!
+! stores <U V>
+IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+ TZFIELD%CMNHNAME = 'UV_FLX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'UV_FLX'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_UV_FLX'
+ TZFIELD%NGRID = 5
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+END IF
+!
+!
+!
+!computation of the source for rho*V due to this flux
+IF (.NOT. LFLAT) THEN
+ PRUS(:,:,:) = PRUS(:,:,:) &
+ - DYF(ZFLX * MXM(MYM(PRHODJ) * PINV_PDYY) ) &
+ + DZF( MYF( MZM(ZFLX)*MXM(PDZY/MZM(PDYY))) &
+ * MXM(PMZM_PRHODJ * PINV_PDZZ) )
+ELSE
+ PRUS(:,:,:) = PRUS(:,:,:) - DYF(ZFLX * MXM(MYM(PRHODJ) * PINV_PDYY) )
+END IF
+!
+!computation of the source for rho*V due to this flux
+IF (.NOT. LFLAT) THEN
+ PRVS(:,:,:) = PRVS(:,:,:) &
+ - DXF(ZFLX * MYM(MXM(PRHODJ) * PINV_PDXX) ) &
+ + DZF( MXF( MZM(ZFLX)*MYM(PDZX/MZM(PDXX))) &
+ * MYM(PMZM_PRHODJ * PINV_PDZZ) )
+ELSE
+ PRVS(:,:,:) = PRVS(:,:,:) - DXF(ZFLX * MYM(MXM(PRHODJ) * PINV_PDXX) )
+END IF
+!
+IF (KSPLT==1) THEN
+ !
+ !Contribution to the dynamic production of TKE:
+ !
+ IF (.NOT. L2D) THEN
+ ZWORK(:,:,:) = - MXF( MYF( ZFLX * &
+ (GY_U_UV_PUM + GX_V_UV_PVM) ) )
+ ELSE
+ ZWORK(:,:,:) = - MXF( MYF( ZFLX * &
+ (GX_V_UV_PVM) ) )
+ ENDIF
+ !
+ ! evaluate the dynamic production at w(IKB+1) in PDP(IKB)
+ !
+ ZWORK(:,:,IKB:IKB) = - &
+ MXF ( MYF( 0.5 * (ZFLX(:,:,IKB+1:IKB+1)+ZFLX(:,:,IKB:IKB)) ) ) &
+ *(MXF ( MYF( &
+ DYM( 0.5 * (PUM(:,:,IKB+1:IKB+1)+PUM(:,:,IKB:IKB)) ) &
+ / MXM( 0.5*(PDYY(:,:,IKB:IKB)+PDYY(:,:,IKB+1:IKB+1)) ) &
+ +DXM( 0.5 * (PVM(:,:,IKB+1:IKB+1)+PVM(:,:,IKB:IKB)) ) &
+ / MYM( 0.5*(PDXX(:,:,IKB:IKB)+PDXX(:,:,IKB+1:IKB+1)) ) &
+ ) ) &
+ -MXF( (PUM(:,:,IKB+1:IKB+1)-PUM(:,:,IKB:IKB)) / &
+ MXM(PDZZ(:,:,IKB+1:IKB+1)) * PDZY(:,:,IKB+1:IKB+1) &
+ ) / MYF(MXM( 0.5*(PDYY(:,:,IKB:IKB)+PDYY(:,:,IKB+1:IKB+1)) ) )&
+ -MYF( (PVM(:,:,IKB+1:IKB+1)-PVM(:,:,IKB:IKB)) / &
+ MYM(PDZZ(:,:,IKB+1:IKB+1)) * PDZX(:,:,IKB+1:IKB+1) &
+ ) / MXF(MYM( 0.5*(PDXX(:,:,IKB:IKB)+PDXX(:,:,IKB+1:IKB+1)) ) )&
+ )
+ !
+ ! dynamic production
+ PDP(:,:,:) = PDP(:,:,:) + ZWORK(:,:,:)
+ !
+END IF
+!
+! Storage in the LES configuration
+!
+IF (LLES_CALL .AND. KSPLT==1) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( MXF(MYF(ZFLX)), X_LES_SUBGRID_UV )
+ CALL LES_MEAN_SUBGRID( MXF(MYF(GY_U_UV(PUM,PDYY,PDZZ,PDZY)*ZFLX)), X_LES_RES_ddxa_U_SBG_UaU , .TRUE.)
+ CALL LES_MEAN_SUBGRID( MXF(MYF(GX_V_UV(PVM,PDXX,PDZZ,PDZX)*ZFLX)), X_LES_RES_ddxa_V_SBG_UaV , .TRUE.)
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+!
+!
+END SUBROUTINE TURB_HOR_UV
+END MODULE MODE_TURB_HOR_UV
diff --git a/src/common/turb/mode_turb_hor_uw.F90 b/src/common/turb/mode_turb_hor_uw.F90
new file mode 100644
index 0000000000000000000000000000000000000000..47005d4a3730f5920f87dc9b813f7f50f6d68985
--- /dev/null
+++ b/src/common/turb/mode_turb_hor_uw.F90
@@ -0,0 +1,248 @@
+!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_HOR_UW
+IMPLICIT NONE
+CONTAINS
+! ################################################################
+ SUBROUTINE TURB_HOR_UW(KSPLT, &
+ OTURB_FLX,KRR, &
+ TPFILE, &
+ PK,PINV_PDXX,PINV_PDZZ,PMZM_PRHODJ, &
+ PDXX,PDZZ,PDZX, &
+ PRHODJ,PTHVREF, &
+ PUM,PWM,PTHLM,PRM,PSVM, &
+ PTKEM,PLM, &
+ PDP, &
+ PRUS,PRWS )
+! ################################################################
+!
+!
+!!**** *TURB_HOR* -routine to compute the source terms in the meso-NH
+!! model equations due to the non-vertical turbulent fluxes.
+!!
+!! PURPOSE
+!! -------
+!!
+!! see TURB_HOR
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!!
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Aug , 1997 (V. Saravane) spliting of TURB_HOR
+!! Nov 27, 1997 (V. Masson) clearing of the routine
+!! Oct 18, 2000 (V. Masson) LES computations + LFLAT switch
+!! Feb 14, 2001 (V. Masson and J. Stein) DZF bug on PRWS
+!! + remove the use of W=0 at the ground
+!! + extrapolation under the ground
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
+!! change of YCOMMENT
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+USE MODD_CONF
+USE MODD_CTURB
+USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_PARAMETERS
+USE MODD_LES
+USE MODD_NSV
+!
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
+!
+USE MODI_GRADIENT_M
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_SHUMAN
+USE MODE_COEFJ, ONLY: COEFJ
+USE MODI_LES_MEAN_SUBGRID
+!
+USE MODI_SECOND_MNH
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declaration of arguments
+!
+!
+!
+INTEGER, INTENT(IN) :: KSPLT ! split process index
+LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
+ ! turbulent fluxes in the syncronous FM-file
+INTEGER, INTENT(IN) :: KRR ! number of moist var.
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PK ! Turbulent diffusion doef.
+ ! PK = PLM * SQRT(PTKEM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDXX ! 1./PDXX
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDZZ ! 1./PDZZ
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PMZM_PRHODJ ! MZM(PRHODJ)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDZZ, PDZX
+ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density * grid volume
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! ref. state VPT
+!
+! Variables at t-1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUM,PWM,PTHLM
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at time t- dt
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRWS
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDP ! TKE production terms
+!
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),SIZE(PWM,3)) &
+ :: ZFLX,ZWORK
+ ! work arrays
+!
+INTEGER :: IKB,IKE,IKU
+ ! Index values for the Beginning and End
+ ! mass points of the domain
+INTEGER :: JSV ! scalar loop counter
+!
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),SIZE(PWM,3)) :: GX_W_UW_PWM
+!
+REAL :: ZTIME1, ZTIME2
+TYPE(TFIELDDATA) :: TZFIELD
+! ---------------------------------------------------------------------------
+!
+!* 1. PRELIMINARY COMPUTATIONS
+! ------------------------
+!
+IKB = 1+JPVEXT
+IKE = SIZE(PWM,3)-JPVEXT
+IKU = SIZE(PWM,3)
+!
+!
+GX_W_UW_PWM = GX_W_UW(PWM,PDXX,PDZZ,PDZX)
+!
+!
+!* 13. < U'W'>
+! -------
+!
+! residual part of < U'W'> depending on dw/dx
+!
+ZFLX(:,:,:) = &
+ - XCMFS * MXM(MZM(PK)) * GX_W_UW_PWM
+!! & to be tested
+!! - (2./3.) * XCMFB * MZM( ZVPTU * MXM( PLM / SQRT(PTKEM) * XG / PTHVREF ) )
+!
+ZFLX(:,:,IKE+1) = 0. ! rigid wall condition => no turbulent flux
+!
+! Nullify the flux at the ground level because it has been fully taken into
+! account in turb_ver and extrapolate the flux under the ground
+ZFLX(:,:,IKB) = 0.
+ZFLX(:,:,IKB-1)=2.*ZFLX(:,:,IKB)- ZFLX(:,:,IKB+1)
+!
+! stores <U W>
+IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+ TZFIELD%CMNHNAME = 'UW_HFLX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'UW_HFLX'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_UW_HFLX'
+ TZFIELD%NGRID = 6
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+END IF
+!
+!
+! compute the source for rho*U due to this residual flux ( the other part is
+! taken into account in TURB_VER)
+PRUS(:,:,:) = PRUS(:,:,:) - DZF( ZFLX* MXM( PMZM_PRHODJ ) / MXM( PDZZ ) )
+!
+!computation of the source for rho*W due to this flux
+IF (.NOT. LFLAT) THEN
+ PRWS(:,:,:) = PRWS(:,:,:) &
+ -DXF( MZM( MXM(PRHODJ) * PINV_PDXX) * ZFLX) &
+ +DZM( PRHODJ * MXF( MZF( ZFLX*PDZX ) * PINV_PDXX ) / MZF(PDZZ) )
+ELSE
+ PRWS(:,:,:) = PRWS(:,:,:) -DXF( MZM( MXM(PRHODJ) * PINV_PDXX) * ZFLX)
+END IF
+!
+IF (KSPLT==1) THEN
+ !
+ !Contribution to the dynamic production of TKE:
+ !
+ ZWORK(:,:,:) =-MZF( MXF( &
+ ZFLX *( GZ_U_UW(PUM,PDZZ) + GX_W_UW_PWM ) ) )
+ !
+ !
+ ! evaluate the dynamic production at w(IKB+1) in PDP(IKB)
+ ZWORK(:,:,IKB:IKB) = - MXF ( &
+ ZFLX(:,:,IKB+1:IKB+1) * &
+ ( (PUM(:,:,IKB+1:IKB+1)-PUM(:,:,IKB:IKB)) / MXM(PDZZ(:,:,IKB+1:IKB+1))&
+ + ( DXM( PWM(:,:,IKB+1:IKB+1) ) &
+ -MXM( (PWM(:,:,IKB+2:IKB+2)-PWM(:,:,IKB+1:IKB+1)) &
+ /(PDZZ(:,:,IKB+2:IKB+2)+PDZZ(:,:,IKB+1:IKB+1)) &
+ +(PWM(:,:,IKB+1:IKB+1)-PWM(:,:,IKB :IKB )) &
+ /(PDZZ(:,:,IKB+1:IKB+1)+PDZZ(:,:,IKB :IKB )) &
+ ) &
+ * PDZX(:,:,IKB+1:IKB+1) &
+ ) / (0.5*(PDXX(:,:,IKB+1:IKB+1)+PDXX(:,:,IKB:IKB))) &
+ ) )
+ !
+ ! dynamic production computation
+ PDP(:,:,:) = PDP(:,:,:) + ZWORK(:,:,:)
+ !
+END IF
+!
+! Storage in the LES configuration (addition to TURB_VER computation)
+!
+IF (LLES_CALL .AND. KSPLT==1) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( MZF(MXF(ZFLX)), X_LES_SUBGRID_WU , .TRUE. )
+ CALL LES_MEAN_SUBGRID( MZF(MXF(GZ_U_UW(PUM,PDZZ)*ZFLX)), X_LES_RES_ddxa_U_SBG_UaU , .TRUE.)
+ CALL LES_MEAN_SUBGRID( MZF(MXF(GX_W_UW_PWM*ZFLX)), X_LES_RES_ddxa_W_SBG_UaW , .TRUE.)
+ CALL LES_MEAN_SUBGRID( MXF(GX_M_U(1,IKU,1,PTHLM,PDXX,PDZZ,PDZX)*MZF(ZFLX)),&
+ X_LES_RES_ddxa_Thl_SBG_UaW , .TRUE.)
+ IF (KRR>=1) THEN
+ CALL LES_MEAN_SUBGRID( MXF(GX_M_U(1,IKU,1,PRM(:,:,:,1),PDXX,PDZZ,PDZX)*MZF(ZFLX)), &
+ X_LES_RES_ddxa_Rt_SBG_UaW , .TRUE.)
+ END IF
+ DO JSV=1,NSV
+ CALL LES_MEAN_SUBGRID( MXF(GX_M_U(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)*MZF(ZFLX)), &
+ X_LES_RES_ddxa_Sv_SBG_UaW(:,:,:,JSV) , .TRUE.)
+ END DO
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+
+!
+END SUBROUTINE TURB_HOR_UW
+END MODULE MODE_TURB_HOR_UW
diff --git a/src/common/turb/mode_turb_hor_vw.F90 b/src/common/turb/mode_turb_hor_vw.F90
new file mode 100644
index 0000000000000000000000000000000000000000..8ede64d515ea715cf2c3f0ee6fc09f7d06909b91
--- /dev/null
+++ b/src/common/turb/mode_turb_hor_vw.F90
@@ -0,0 +1,259 @@
+!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_HOR_VW
+IMPLICIT NONE
+CONTAINS
+ SUBROUTINE TURB_HOR_VW(KSPLT, &
+ OTURB_FLX,KRR, &
+ TPFILE, &
+ PK,PINV_PDYY,PINV_PDZZ,PMZM_PRHODJ, &
+ PDYY,PDZZ,PDZY, &
+ PRHODJ,PTHVREF, &
+ PVM,PWM,PTHLM,PRM,PSVM, &
+ PTKEM,PLM, &
+ PDP, &
+ PRVS,PRWS )
+! ################################################################
+!
+!
+!!**** *TURB_HOR* -routine to compute the source terms in the meso-NH
+!! model equations due to the non-vertical turbulent fluxes.
+!!
+!! PURPOSE
+!! -------
+!!
+!! see TURB_HOR
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!!
+!! Joan Cuxart * INM and Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Aug , 1997 (V. Saravane) spliting of TURB_HOR
+!! Nov 27, 1997 (V. Masson) clearing of the routine
+!! Oct 18, 2000 (V. Masson) LES computations + LFLAT switch
+!! Feb 14, 2001 (V. Masson and J. Stein) DZF bug on PRWS
+!! + remove the use of W=0 at the ground
+!! + extrapolataion under the ground
+!! Nov 06, 2002 (V. Masson) LES budgets
+!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
+!! change of YCOMMENT
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+USE MODD_CONF
+USE MODD_CTURB
+USE MODD_FIELD, ONLY: TFIELDDATA, TYPEREAL
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_PARAMETERS
+USE MODD_LES
+USE MODD_NSV
+!
+USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
+!
+USE MODI_GRADIENT_M
+USE MODI_GRADIENT_U
+USE MODI_GRADIENT_V
+USE MODI_GRADIENT_W
+USE MODI_SHUMAN
+USE MODE_COEFJ, ONLY: COEFJ
+USE MODI_LES_MEAN_SUBGRID
+!
+USE MODI_SECOND_MNH
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declaration of arguments
+!
+!
+!
+INTEGER, INTENT(IN) :: KSPLT ! split process index
+LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
+ ! turbulent fluxes in the syncronous FM-file
+INTEGER, INTENT(IN) :: KRR ! number of moist var.
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PK ! Turbulent diffusion doef.
+ ! PK = PLM * SQRT(PTKEM)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDYY ! 1./PDYY
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDZZ ! 1./PDZZ
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PMZM_PRHODJ ! MZM(PRHODJ)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY, PDZZ, PDZY
+ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density * grid volume
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! ref. state VPT
+!
+! Variables at t-1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PVM,PWM,PTHLM
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at time t- dt
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRVS, PRWS ! var. at t+1 -split-
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDP ! TKE production terms
+!
+!
+!
+!* 0.2 declaration of local variables
+!
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),SIZE(PWM,3)) &
+ :: ZFLX,ZWORK
+ ! work arrays
+!
+!! REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),SIZE(PWM,3)) :: ZVPTV
+INTEGER :: IKB,IKE,IKU
+ ! Index values for the Beginning and End
+ ! mass points of the domain
+INTEGER :: JSV ! scalar loop counter
+!
+REAL, DIMENSION(SIZE(PWM,1),SIZE(PWM,2),SIZE(PWM,3)) :: GY_W_VW_PWM
+!
+REAL :: ZTIME1, ZTIME2
+TYPE(TFIELDDATA) :: TZFIELD
+! ---------------------------------------------------------------------------
+!
+!* 1. PRELIMINARY COMPUTATIONS
+! ------------------------
+!
+IKB = 1+JPVEXT
+IKE = SIZE(PWM,3)-JPVEXT
+IKU = SIZE(PWM,3)
+!
+!
+IF (.NOT. L2D) GY_W_VW_PWM = GY_W_VW(PWM,PDYY,PDZZ,PDZY)
+!
+!
+!* 14. < V'W'>
+! -------
+!
+! residual part of < V'W'> depending on dw/dy
+!
+IF (.NOT. L2D) THEN
+ ZFLX(:,:,:) = &
+ - XCMFS * MYM(MZM(PK)) * GY_W_VW_PWM
+ !! & to be tested
+ !! - (2./3.) * XCMFB * MZM( ZVPTV * MYM( PLM / SQRT(PTKEM) * XG / PTHVREF ) )
+ELSE
+ ZFLX(:,:,:) = 0.
+ !! & to be tested
+ !! - (2./3.) * XCMFB * MZM( ZVPTV * MYM( PLM / SQRT(PTKEM) * XG / PTHVREF ) )
+END IF
+!
+ZFLX(:,:,IKE+1) = 0. ! rigid wall condition => no turbulent flux
+!
+!
+! Nullify the flux at the ground level because it has been fully taken into
+! account in turb_ver and extrapolate the flux under the ground
+ZFLX(:,:,IKB) = 0.
+ZFLX(:,:,IKB-1)= 2.*ZFLX(:,:,IKB) - ZFLX(:,:,IKB+1)
+!
+! stores <V W>
+IF ( TPFILE%LOPENED .AND. OTURB_FLX ) THEN
+ TZFIELD%CMNHNAME = 'VW_HFLX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'VW_HFLX'
+ TZFIELD%CUNITS = 'm2 s-2'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_VW_HFLX'
+ TZFIELD%NGRID = 7
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZFLX)
+END IF
+!
+! compute the source for rho*V due to this residual flux ( the other part is
+! taken into account in TURB_VER)
+IF (.NOT. L2D) &
+PRVS(:,:,:) = PRVS(:,:,:) - DZF( ZFLX* MYM( PMZM_PRHODJ ) / MYM ( PDZZ ) )
+!
+!computation of the source for rho*W due to this flux
+IF (.NOT. L2D) THEN
+ IF (.NOT. LFLAT) THEN
+ PRWS(:,:,:) = PRWS(:,:,:) &
+ -DYF( MZM( MYM(PRHODJ) * PINV_PDYY) * ZFLX) &
+ +DZM( PRHODJ * MYF( MZF( ZFLX*PDZY ) * PINV_PDYY ) / MZF(PDZZ) )
+ ELSE
+ PRWS(:,:,:) = PRWS(:,:,:) - DYF( MZM( MYM(PRHODJ) * PINV_PDYY) * ZFLX)
+ END IF
+END IF
+!
+IF (KSPLT==1) THEN
+ !
+ !Contribution to the dynamic production of TKE:
+ !
+ IF (.NOT. L2D) THEN
+ ZWORK(:,:,:) =-MZF( MYF( ZFLX *( GZ_V_VW(PVM,PDZZ) + GY_W_VW_PWM ) ) )
+ !
+ !
+ ! evaluate the dynamic production at w(IKB+1) in PDP(IKB)
+ ZWORK(:,:,IKB:IKB) = - MYF ( &
+ ZFLX(:,:,IKB+1:IKB+1) * &
+ ( (PVM(:,:,IKB+1:IKB+1)-PVM(:,:,IKB:IKB)) / MYM(PDZZ(:,:,IKB+1:IKB+1)) &
+ + ( DYM( PWM(:,:,IKB+1:IKB+1) ) &
+ -MYM( (PWM(:,:,IKB+2:IKB+2)-PWM(:,:,IKB+1:IKB+1)) &
+ /(PDZZ(:,:,IKB+2:IKB+2)+PDZZ(:,:,IKB+1:IKB+1)) &
+ +(PWM(:,:,IKB+1:IKB+1)-PWM(:,:,IKB :IKB )) &
+ /(PDZZ(:,:,IKB+1:IKB+1)+PDZZ(:,:,IKB :IKB )) &
+ ) * PDZY(:,:,IKB+1:IKB+1) &
+ ) / (0.5*(PDYY(:,:,IKB+1:IKB+1)+PDYY(:,:,IKB:IKB))) &
+ ) )
+ ENDIF
+ !
+ ! dynamic production computation
+ IF (.NOT. L2D) &
+ PDP(:,:,:) = PDP(:,:,:) + ZWORK(:,:,:)
+ !
+END IF
+!
+! Storage in the LES configuration (addition to TURB_VER computation)
+!
+IF (LLES_CALL .AND. KSPLT==1) THEN
+ CALL SECOND_MNH(ZTIME1)
+ CALL LES_MEAN_SUBGRID( MZF(MYF(ZFLX)), X_LES_SUBGRID_WV , .TRUE. )
+ CALL LES_MEAN_SUBGRID( MZF(MYF(GZ_V_VW(PVM,PDZZ)*ZFLX)),&
+ X_LES_RES_ddxa_V_SBG_UaV , .TRUE.)
+ CALL LES_MEAN_SUBGRID( MZF(MYF(GY_W_VW(PWM,PDYY,PDZZ,PDZY)*ZFLX)),&
+ X_LES_RES_ddxa_W_SBG_UaW , .TRUE.)
+ CALL LES_MEAN_SUBGRID( MXF(GY_M_V(1,IKU,1,PTHLM,PDYY,PDZZ,PDZY)*MZF(ZFLX)),&
+ X_LES_RES_ddxa_Thl_SBG_UaW , .TRUE.)
+ IF (KRR>=1) THEN
+ CALL LES_MEAN_SUBGRID( MXF(GY_M_V(1,IKU,1,PRM(:,:,:,1),PDYY,PDZZ,PDZY)*MZF(ZFLX)), &
+ X_LES_RES_ddxa_Rt_SBG_UaW , .TRUE.)
+ END IF
+ DO JSV=1,NSV
+ CALL LES_MEAN_SUBGRID( MXF(GY_M_V(1,IKU,1,PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)*MZF(ZFLX)), &
+ X_LES_RES_ddxa_Sv_SBG_UaW(:,:,:,JSV), .TRUE.)
+ END DO
+ CALL SECOND_MNH(ZTIME2)
+ XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
+END IF
+!
+!
+!
+END SUBROUTINE TURB_HOR_VW
+END MODULE MODE_TURB_HOR_VW
diff --git a/src/common/turb/turb.F90 b/src/common/turb/turb.F90
index 989437065c9e417ea4f2ccbaaaa3ee7d3c826fb0..2e97d63c673ebe1703d8079c6907c4b0903eb740 100644
--- a/src/common/turb/turb.F90
+++ b/src/common/turb/turb.F90
@@ -242,7 +242,7 @@ USE MODE_BL89, ONLY: BL89
USE MODE_TURB_VER, ONLY : TURB_VER
!!MODIF AROME
!USE MODI_ROTATE_WIND
-!USE MODE_TURB_HOR_SPLT, ONLY: TURB_HOR_SPLT
+USE MODE_TURB_HOR_SPLT, ONLY: TURB_HOR_SPLT
USE MODE_TKE_EPS_SOURCES, ONLY: TKE_EPS_SOURCES
USE MODI_SHUMAN, ONLY : MZF, MXF, MYF
USE MODI_GRADIENT_M
@@ -985,24 +985,24 @@ IF( HTURBDIM == '3DIM' ) THEN
CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'HTURB', PRSVS(:, :, :, JSV) )
END DO
END IF
-!#ifdef REPRO48 !à supprimer une fois le précédent ifdef REPRO48 validé
-!#else
-! CALL TURB_HOR_SPLT(KSPLIT, KRR, KRRL, KRRI, PTSTEP, &
-! HLBCX,HLBCY,OTURB_FLX,OSUBG_COND, &
-! TPFILE, &
-! PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, &
-! PDIRCOSXW,PDIRCOSYW,PDIRCOSZW, &
-! PCOSSLOPE,PSINSLOPE, &
-! PRHODJ,PTHVREF, &
-! PSFTH,PSFRV,PSFSV, &
-! ZCDUEFF,ZTAU11M,ZTAU12M,ZTAU22M,ZTAU33M, &
-! PUT,PVT,PWT,ZUSLOPE,ZVSLOPE,PTHLT,PRT,PSVT, &
-! PTKET,PLEM,ZLEPS, &
-! ZLOCPEXNM,ZATHETA,ZAMOIST,PSRCT,ZFRAC_ICE, &
-! PDYP,PTHP,PSIGS, &
-! ZTRH, &
-! PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS )
-!#endif
+#ifdef REPRO48 !à supprimer une fois le précédent ifdef REPRO48 validé
+#else
+ CALL TURB_HOR_SPLT(KSPLIT, KRR, KRRL, KRRI, PTSTEP, &
+ HLBCX,HLBCY,OTURB_FLX,OSUBG_COND, &
+ TPFILE, &
+ PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, &
+ PDIRCOSXW,PDIRCOSYW,PDIRCOSZW, &
+ PCOSSLOPE,PSINSLOPE, &
+ PRHODJ,PTHVREF, &
+ PSFTH,PSFRV,PSFSV, &
+ ZCDUEFF,ZTAU11M,ZTAU12M,ZTAU22M,ZTAU33M, &
+ PUT,PVT,PWT,ZUSLOPE,ZVSLOPE,PTHLT,PRT,PSVT, &
+ PTKET,PLEM,ZLEPS, &
+ ZLOCPEXNM,ZATHETA,ZAMOIST,PSRCT,ZFRAC_ICE, &
+ PDYP,PTHP,PSIGS, &
+ ZTRH, &
+ PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS )
+#endif
IF( LBUDGET_U ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_U), 'HTURB', PRUS(:, :, :) )
IF( LBUDGET_V ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_V), 'HTURB', PRVS(:, :, :) )
IF( LBUDGET_W ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_W), 'HTURB', PRWS(:, :, :) )