From 91915bdeb61d0b0dfea13b8e14bdb1be7b3e535b Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Tue, 8 Jun 2021 11:57:05 +0200
Subject: [PATCH] Didier R. 08/06/2021: add Leonard terms (Moeng
parametrization) instead of Reynolds terms applied to vertical fluxes of r_np
and Thl for implicit version of turbulence scheme
---
src/MNH/default_desfmn.f90 | 8 ++-
src/MNH/modd_turbn.f90 | 18 ++++++
src/MNH/modn_turbn.f90 | 26 ++++++++-
src/MNH/turb_ver_thermo_flux.f90 | 97 ++++++++++++++++++++++++++++++--
4 files changed, 141 insertions(+), 8 deletions(-)
diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index 960f34cbf..f39f69423 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -212,7 +212,7 @@ END MODULE MODI_DEFAULT_DESFM_n
!! 02/2021 (T.Nagel) add turbulence recycling defaults parameters
! P-A Joulin 21/05/2021: add Wind turbines
! S. Riette 21/05/2021: add options to PDF subgrid scheme
-!
+! 05/2021 D. Ricard add the contribution of Leonard terms in the turbulence scheme
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -533,6 +533,12 @@ VSIGQSAT = 0.02
CCONDENS='CB02'
CLAMBDA3='CB'
CSUBG_MF_PDF='TRIANGLE'
+LHGRAD =.FALSE.
+XCOEFHGRADTHL = 1.0
+XCOEFHGRADRM = 1.0
+XALTHGRAD = 2000.0
+XCLDTHOLD = -1.0
+
!-------------------------------------------------------------------------------
!
!* 10b. SET DEFAULT VALUES FOR MODD_DRAGTREE :
diff --git a/src/MNH/modd_turbn.f90 b/src/MNH/modd_turbn.f90
index a108cce8f..ac0b49ae3 100644
--- a/src/MNH/modd_turbn.f90
+++ b/src/MNH/modd_turbn.f90
@@ -39,6 +39,7 @@
!! May 2006 Remove KEPS
!! C.Lac Nov 2014 add terms of TKE production for LES diag
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! D. Ricard May 2021 add the switches for Leonard terms
!!
!-------------------------------------------------------------------------------
!
@@ -93,6 +94,13 @@ TYPE TURB_t
REAL, DIMENSION(:,:,:), POINTER :: XTR=>NULL() ! Transport production of Kinetic energy
REAL, DIMENSION(:,:,:), POINTER :: XDISS=>NULL() ! Dissipation of Kinetic energy
REAL, DIMENSION(:,:,:), POINTER :: XLEM=>NULL() ! Mixing length
+ LOGICAL :: LHGRAD ! logical switch for the computation of the Leornard Terms
+ REAL :: XCOEFHGRADTHL ! coeff applied to thl contribution
+ REAL :: XCOEFHGRADRM ! coeff applied to mixing ratio contribution
+ REAL :: XALTHGRAD ! altitude from which to apply the Leonard terms
+ REAL :: XCLDTHOLD ! cloud threshold to apply the Leonard terms
+ ! negative value : applied everywhere
+ ! 0.000001 applied only inside the clouds ri+rc > 10**-6 kg/kg
!
END TYPE TURB_t
@@ -125,6 +133,11 @@ REAL, DIMENSION(:,:,:), POINTER :: XTHP=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XTR=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XDISS=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XLEM=>NULL()
+LOGICAL, POINTER :: LHGRAD=>NULL()
+REAL, POINTER :: XCOEFHGRADTHL=>NULL()
+REAL, POINTER :: XCOEFHGRADRM=>NULL()
+REAL, POINTER :: XALTHGRAD=>NULL()
+REAL, POINTER :: XCLDTHOLD=>NULL()
CONTAINS
@@ -170,6 +183,11 @@ XTHP=>TURB_MODEL(KTO)%XTHP
XTR=>TURB_MODEL(KTO)%XTR
XDISS=>TURB_MODEL(KTO)%XDISS
XLEM=>TURB_MODEL(KTO)%XLEM
+LHGRAD=>TURB_MODEL(KTO)%LHGRAD
+XCOEFHGRADTHL=>TURB_MODEL(KTO)%XCOEFHGRADTHL
+XCOEFHGRADRM=>TURB_MODEL(KTO)%XCOEFHGRADRM
+XALTHGRAD=>TURB_MODEL(KTO)%XALTHGRAD
+XCLDTHOLD=>TURB_MODEL(KTO)%XCLDTHOLD
END SUBROUTINE TURB_GOTO_MODEL
diff --git a/src/MNH/modn_turbn.f90 b/src/MNH/modn_turbn.f90
index 96a4b347c..952412dbe 100644
--- a/src/MNH/modn_turbn.f90
+++ b/src/MNH/modn_turbn.f90
@@ -47,6 +47,7 @@
!! P. Bechtold Feb 11, 2002 add switch for Sigma_s computation
!! P. Jabouille Apr 4, 2002 add switch for Sigma_s convection
!! V. Masson Nov 13 2002 add switch for SBL lengths
+!! D. Ricard May, 2021 add switch for Leonard Terms
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -71,7 +72,12 @@ USE MODD_TURB_n, ONLY: &
CSUBG_AUCV_RI_n => CSUBG_AUCV_RI, &
CCONDENS_n => CCONDENS, &
CLAMBDA3_n => CLAMBDA3, &
- CSUBG_MF_PDF_n => CSUBG_MF_PDF
+ CSUBG_MF_PDF_n => CSUBG_MF_PDF, &
+ LHGRAD_n => LHGRAD, &
+ XCOEFHGRADTHL_n => XCOEFHGRADTHL, &
+ XCOEFHGRADRM_n => XCOEFHGRADRM, &
+ XALTHGRAD_n => XALTHGRAD, &
+ XCLDTHOLD_n => XCLDTHOLD
!
IMPLICIT NONE
!
@@ -94,11 +100,17 @@ CHARACTER (LEN=80),SAVE :: CCONDENS
CHARACTER (LEN=4),SAVE :: CLAMBDA3
CHARACTER (LEN=80),SAVE :: CSUBG_MF_PDF
REAL,SAVE :: VSIGQSAT
+LOGICAL,SAVE :: LHGRAD
+REAL,SAVE :: XCOEFHGRADTHL
+REAL,SAVE :: XCOEFHGRADRM
+REAL,SAVE :: XALTHGRAD
+REAL,SAVE :: XCLDTHOLD
!
NAMELIST/NAM_TURBn/XIMPL,CTURBLEN,CTURBDIM,LTURB_FLX,LTURB_DIAG, &
LSUBG_COND,LSIGMAS,LSIG_CONV,LRMC01,CTOM,CSUBG_AUCV,&
XKEMIN,VSIGQSAT,XCEDIS,XCADAP,CSUBG_AUCV_RI,CCONDENS,&
- CLAMBDA3,CSUBG_MF_PDF
+ CLAMBDA3,CSUBG_MF_PDF,LHGRAD,XCOEFHGRADTHL, XCOEFHGRADRM, &
+ XALTHGRAD, XCLDTHOLD
!
CONTAINS
@@ -123,6 +135,11 @@ SUBROUTINE INIT_NAM_TURBn
CCONDENS = CCONDENS_n
CLAMBDA3 = CLAMBDA3_n
CSUBG_MF_PDF = CSUBG_MF_PDF_n
+ LHGRAD = LHGRAD_n
+ XCOEFHGRADTHL = XCOEFHGRADTHL_n
+ XCOEFHGRADRM = XCOEFHGRADRM_n
+ XALTHGRAD = XALTHGRAD_n
+ XCLDTHOLD = XCLDTHOLD_n
END SUBROUTINE INIT_NAM_TURBn
SUBROUTINE UPDATE_NAM_TURBn
@@ -145,6 +162,11 @@ SUBROUTINE UPDATE_NAM_TURBn
CCONDENS_n = CCONDENS
CLAMBDA3_n = CLAMBDA3
CSUBG_MF_PDF_n = CSUBG_MF_PDF
+ LHGRAD_n = LHGRAD
+ XCOEFHGRADTHL_n = XCOEFHGRADTHL
+ XCOEFHGRADRM_n = XCOEFHGRADRM
+ XALTHGRAD_n = XALTHGRAD
+ XCLDTHOLD_n = XCLDTHOLD
END SUBROUTINE UPDATE_NAM_TURBn
END MODULE MODN_TURB_n
diff --git a/src/MNH/turb_ver_thermo_flux.f90 b/src/MNH/turb_ver_thermo_flux.f90
index 80825e315..22fdeae7c 100644
--- a/src/MNH/turb_ver_thermo_flux.f90
+++ b/src/MNH/turb_ver_thermo_flux.f90
@@ -323,6 +323,11 @@ END MODULE MODI_TURB_VER_THERMO_FLUX
!! 2012-02 (Y. Seity) add possibility to run with reversed
!! vertical levels
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! 2018 (D. Ricard) last version of HGRAD turbulence scheme
+!! Leronard terms instead of Reynolds terms
+!! applied to vertical fluxes of r_np and Thl
+!! for implicit version of turbulence scheme
+!! corrections and cleaning
!!--------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -331,8 +336,11 @@ END MODULE MODI_TURB_VER_THERMO_FLUX
USE MODD_CST
USE MODD_CTURB
use modd_field, only: tfielddata, TYPEREAL
+USE MODD_GRID_n, ONLY: XZS, XXHAT
USE MODD_IO, ONLY: TFILEDATA
+USE MODD_METRICS_n
USE MODD_PARAMETERS
+USE MODD_TURB_n, ONLY: LHGRAD, XCOEFHGRADTHL, XCOEFHGRADRM, XALTHGRAD, XCLDTHOLD
USE MODD_CONF
USE MODD_LES
!
@@ -340,6 +348,9 @@ USE MODI_GRADIENT_U
USE MODI_GRADIENT_V
USE MODI_GRADIENT_W
USE MODI_GRADIENT_M
+USE MODI_GRADIENT_UV
+USE MODI_GRADIENT_UW
+USE MODI_GRADIENT_VW
USE MODI_SHUMAN
USE MODI_TRIDIAG
USE MODI_LES_MEAN_SUBGRID
@@ -452,13 +463,21 @@ REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) :: &
ZF, & ! Flux in dTh/dt =-dF/dz (evaluated at t-1)(or rt instead of Th)
ZDFDDTDZ, & ! dF/d(dTh/dz)
ZDFDDRDZ, & ! dF/d(dr/dz)
- Z3RDMOMENT ! 3 order term in flux or variance equation
+ Z3RDMOMENT,& ! 3 order term in flux or variance equation
+ ZF_NEW, &
+ ZRWTHL, &
+ ZRWRNP, &
+ ZCLD_THOLD
+!
+REAL,DIMENSION(SIZE(XZS,1),SIZE(XZS,2),KKU) :: ZALT
+!
INTEGER :: IKB,IKE ! I index values for the Beginning and End
! mass points of the domain in the 3 direct.
INTEGER :: IKT ! array size in k direction
INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
!
REAL :: ZTIME1, ZTIME2
+REAL :: ZDELTAX
!
INTEGER :: JK
LOGICAL :: GUSERV ! flag to use water
@@ -486,6 +505,18 @@ GUSERV = (KRR/=0)
!
ZKEFF(:,:,:) = MZM( PLM(:,:,:) * SQRT(PTKEM(:,:,:)) )
!
+! define a cloud mask with ri and rc (used after with a threshold) for Leonard terms
+!
+IF(LHGRAD) THEN
+ IF ( KRRL >= 1 ) THEN
+ IF ( KRRI >= 1 ) THEN
+ ZCLD_THOLD(:,:,:) = PRM(:,:,:,2) + PRM(:,:,:,4)
+ ELSE
+ ZCLD_THOLD(:,:,:) = PRM(:,:,:,2)
+ END IF
+ END IF
+END IF
+!
! Flags for 3rd order quantities
!
GFTH2 = .FALSE.
@@ -514,6 +545,16 @@ END IF
ZF (:,:,:) = -XCSHF*PPHI3*ZKEFF*DZM(PTHLM)/PDZZ
ZDFDDTDZ(:,:,:) = -XCSHF*ZKEFF*D_PHI3DTDZ_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV)
+!
+IF (LHGRAD) THEN
+ ! Compute the Leonard terms for thl
+ ZDELTAX= XXHAT(3) - XXHAT(2)
+ ZF_NEW (:,:,:)= XCOEFHGRADTHL*ZDELTAX*ZDELTAX/12.0*( &
+ MXF(GX_W_UW(PWM(:,:,:), XDXX, XDZZ, XDZX))&
+ *MZM(GX_M_M(PTHLM(:,:,:),XDXX,XDZZ,XDZX)) &
+ + MYF(GY_W_VW(PWM(:,:,:), XDYY,XDZZ,XDZY)) &
+ *MZM(GY_M_M(PTHLM(:,:,:),XDYY,XDZZ,XDZY)) )
+END IF
!
! Effect of 3rd order terms in temperature flux (at flux point)
!
@@ -581,8 +622,19 @@ CALL TRIDIAG_THERMO(KKA,KKU,KKL,PTHLM,ZF,ZDFDDTDZ,PTSTEP,PIMPL,PDZZ,&
PRHODJ,PTHLP)
!
! Compute the equivalent tendency for the conservative potential temperature
-PRTHLS(:,:,:)= PRTHLS(:,:,:) + &
- PRHODJ(:,:,:)*(PTHLP(:,:,:)-PTHLM(:,:,:))/PTSTEP
+!
+ZRWTHL(:,:,:)= PRHODJ(:,:,:)*(PTHLP(:,:,:)-PTHLM(:,:,:))/PTSTEP
+! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
+IF (LHGRAD) THEN
+ DO JK=1,KKU
+ ZALT(:,:,JK) = PZZ(:,:,JK)-XZS(:,:)
+ END DO
+ WHERE ( (ZCLD_THOLD(:,:,:) >= XCLDTHOLD) .AND. ( ZALT(:,:,:) >= XALTHGRAD) )
+ ZRWTHL(:,:,:) = -GZ_W_M(MZM(PRHODJ(:,:,:))*ZF_NEW(:,:,:),XDZZ)
+ END WHERE
+END IF
+!
+PRTHLS(:,:,:)= PRTHLS(:,:,:) + ZRWTHL(:,:,:)
!
!* 2.2 Partial Thermal Production
!
@@ -590,6 +642,12 @@ PRTHLS(:,:,:)= PRTHLS(:,:,:) + &
!
ZFLXZ(:,:,:) = ZF &
+ PIMPL * ZDFDDTDZ * DZM(PTHLP - PTHLM) / PDZZ
+! replace the flux by the Leonard terms
+IF (LHGRAD) THEN
+ WHERE ( (ZCLD_THOLD(:,:,:) >= XCLDTHOLD) .AND. ( ZALT(:,:,:) >= XALTHGRAD) )
+ ZFLXZ(:,:,:) = ZF_NEW(:,:,:)
+ END WHERE
+END IF
!
ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
!
@@ -693,6 +751,16 @@ IF (KRR /= 0) THEN
ZF (:,:,:) = -XCSHF*PPSI3*ZKEFF*DZM(PRM(:,:,:,1))/PDZZ
ZDFDDRDZ(:,:,:) = -XCSHF*ZKEFF*D_PSI3DRDZ_O_DDRDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV)
!
+ ! Compute Leonard Terms for Cloud mixing ratio
+ IF (LHGRAD) THEN
+ ZDELTAX= XXHAT(3) - XXHAT(2)
+ ZF_NEW (:,:,:)= XCOEFHGRADRM*ZDELTAX*ZDELTAX/12.0*( &
+ MXF(GX_W_UW(PWM(:,:,:), XDXX, XDZZ, XDZX)) &
+ *MZM(GX_M_M(PRM(:,:,:,1),XDXX,XDZZ,XDZX)) &
+ +MYF(GY_W_VW(PWM(:,:,:), XDYY,XDZZ,XDZY)) &
+ *MZM(GY_M_M(PRM(:,:,:,1),XDYY,XDZZ,XDZY)) )
+ END IF
+ !
! Effect of 3rd order terms in temperature flux (at flux point)
!
! d(w'2r')/dz
@@ -759,8 +827,20 @@ IF (KRR /= 0) THEN
PDZZ,PRHODJ,PRP)
!
! Compute the equivalent tendency for the conservative mixing ratio
- PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRHODJ(:,:,:) * &
- (PRP(:,:,:)-PRM(:,:,:,1))/PTSTEP
+ !
+ ZRWRNP (:,:,:) = PRHODJ(:,:,:)*(PRP(:,:,:)-PRM(:,:,:,1))/PTSTEP
+ !
+ ! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
+ IF (LHGRAD) THEN
+ DO JK=1,KKU
+ ZALT(:,:,JK) = PZZ(:,:,JK)-XZS(:,:)
+ END DO
+ WHERE ( (ZCLD_THOLD(:,:,:) >= XCLDTHOLD ) .AND. ( ZALT(:,:,:) >= XALTHGRAD ) )
+ ZRWRNP (:,:,:) = -GZ_W_M(MZM(PRHODJ(:,:,:))*ZF_NEW(:,:,:),XDZZ)
+ END WHERE
+ END IF
+ !
+ PRRS(:,:,:,1) = PRRS(:,:,:,1) + ZRWRNP (:,:,:)
!
!* 3.2 Complete thermal production
!
@@ -769,6 +849,13 @@ IF (KRR /= 0) THEN
ZFLXZ(:,:,:) = ZF &
+ PIMPL * ZDFDDRDZ * DZM(PRP - PRM(:,:,:,1)) / PDZZ
!
+ ! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
+ IF (LHGRAD) THEN
+ WHERE ( (ZCLD_THOLD(:,:,:) >= XCLDTHOLD ) .AND. ( ZALT(:,:,:) >= XALTHGRAD ) )
+ ZFLXZ(:,:,:) = ZF_NEW(:,:,:)
+ END WHERE
+ END IF
+ !
ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
!
DO JK=IKTB+1,IKTE-1
--
GitLab