diff --git a/src/arome/ext/aro_turb_mnh.h b/src/arome/ext/aro_turb_mnh.h
index 87d261e866fd200d4e6b43fc953da7c76ac0b63b..bf20e5b9bdaaec12ecdce19101f7c1559dedc09b 100644
--- a/src/arome/ext/aro_turb_mnh.h
+++ b/src/arome/ext/aro_turb_mnh.h
@@ -69,7 +69,7 @@ REAL(KIND=JPRB), DIMENSION(KLON,1,KLEV+2), INTENT(OUT) :: PDRTHLS_TURB
REAL(KIND=JPRB), DIMENSION(KLON,1,KLEV+2), INTENT(OUT) :: PDRRTS_TURB
REAL(KIND=JPRB), DIMENSION(KLON,1,KLEV,KSV), INTENT(OUT) :: PDRSVS_TURB
REAL(KIND=JPRB), DIMENSION(KLON,1,KLEV+2), INTENT(INOUT) :: PFLXZTHVMF
-LOGICAL , INTENT(INOUT) :: OSUBG_COND
+LOGICAL , INTENT(IN) :: OSUBG_COND
REAL(KIND=JPRB), DIMENSION(KLON,1,KLEV+2), INTENT(OUT) :: PDP
REAL(KIND=JPRB), DIMENSION(KLON,1,KLEV+2), INTENT(OUT) :: PTP
REAL(KIND=JPRB), DIMENSION(KLON,1,KLEV+2), INTENT(OUT) :: PTPMF
diff --git a/src/arome/ext/aroini_turb.F90 b/src/arome/ext/aroini_turb.F90
index 6a9edd9d97e6a27dc15a17f06d2be8d18a759f3f..c1ab16006739d7ac4bcb4a2c8d5dd1937d2ab356 100644
--- a/src/arome/ext/aroini_turb.F90
+++ b/src/arome/ext/aroini_turb.F90
@@ -1,5 +1,5 @@
! ######spl
-SUBROUTINE AROINI_TURB(PLINI,OHARATU)
+SUBROUTINE AROINI_TURB(PLINI,OHARATU,OSTATNW)
USE PARKIND1, ONLY : JPRB
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
!**** *INI_TURB* - Initialize common meso_NH MODD_ used in Turbulence scheme
@@ -42,7 +42,7 @@ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
! ------------------------------------------------------------------
USE MODD_LES, ONLY : LLES, LLES_CALL
-USE MODD_CTURB, ONLY : XLINI, LHARAT
+USE MODD_CTURB, ONLY : XLINI, LHARAT,LSTATNW
USE MODD_TURB_n, ONLY : CTURBLEN
USE MODI_INI_CTURB
@@ -54,6 +54,7 @@ IMPLICIT NONE
!
REAL, INTENT(IN) :: PLINI ! minimum bl89 mixing length
LOGICAL,INTENT(IN) :: OHARATU ! switch HARATU
+LOGICAL,INTENT(IN) :: OSTATNW ! switch LSTATNW
!
! ------------------------------------------------------------------
@@ -66,6 +67,7 @@ CALL INI_CTURB
! 1bis. Modification of MODD_CTURB values
XLINI=PLINI
LHARAT=OHARATU
+LSTATNW=OSTATNW
! 2. Set implicit default values for MODD_LES
diff --git a/src/arome/ext/aroini_turb.h b/src/arome/ext/aroini_turb.h
new file mode 100644
index 0000000000000000000000000000000000000000..eeb0a5ab3191974f4ce769e92768c6a5cbe97e63
--- /dev/null
+++ b/src/arome/ext/aroini_turb.h
@@ -0,0 +1,8 @@
+INTERFACE
+SUBROUTINE AROINI_TURB(PLINI,OHARATU,OSTATNW)
+USE PARKIND1 ,ONLY : JPRB
+REAL(KIND=JPRB),INTENT(IN)::PLINI
+LOGICAL, INTENT(IN) ::OHARATU
+LOGICAL, INTENT(IN) ::OSTATNW
+END SUBROUTINE AROINI_TURB
+END INTERFACE
diff --git a/src/arome/ext/suphmpa.F90 b/src/arome/ext/suphmpa.F90
index c5ce5a44ed1c82a1a307712563733980545a4ea9..d5817d3a4b916a46b4d026e2d0b774eeea5e4350 100644
--- a/src/arome/ext/suphmpa.F90
+++ b/src/arome/ext/suphmpa.F90
@@ -160,7 +160,8 @@ CALL AROINI_BUDGET(LAROBU_ENABLE)
! 4. Initialisation of Turbulence scheme
-CALL AROINI_TURB(XLINI,LHARATU)
+!SR phasing step, .FALSE. must be replaced by the value read in namelist
+CALL AROINI_TURB(XLINI,LHARATU,.FALSE.)
! 5. Initialisation of Mass Flux Shallow convection scheme
diff --git a/src/arome/turb/ini_cturb.F90 b/src/arome/turb/ini_cturb.F90
index 1847d0fdaa16277c181230c08f4b43093856b0a7..8880cce0bdda185c5af3f11e65656d457a1cb2a0 100644
--- a/src/arome/turb/ini_cturb.F90
+++ b/src/arome/turb/ini_cturb.F90
@@ -65,7 +65,6 @@ IF (LHOOK) CALL DR_HOOK('INI_CTURB',0,ZHOOK_HANDLE)
! 1.1 Constant for dissipation of Tke
!
!
-LHARAT=.FALSE.
!
!XCED = 0.70
XCED = 0.85
@@ -118,7 +117,12 @@ XCTD = 1.2
!
! 1.7 Constant for temperature and vapor pressure-correlations
!
-XCTP = 4.65
+!wc in STATNW consistent use of Redelsperger-Sommeria for (co)variances
+IF (LSTATNW) THEN
+ XCTP = 4.0
+ ELSE
+ XCTP = 4.65
+ENDIF
! Redelsperger-Sommeria (1981) = 4.
! Schmidt-Schumann (1989) = 3.25
! Cheng-Canuto-Howard (2002) = 4.65
diff --git a/src/common/turb/modd_cturb.F90 b/src/common/turb/modd_cturb.F90
index a446914f779a1f5aa6d72c4c24b7b9a2960261e5..75d4a5948b9a10399fd04e986cef957f21050856 100644
--- a/src/common/turb/modd_cturb.F90
+++ b/src/common/turb/modd_cturb.F90
@@ -82,5 +82,6 @@ REAL,SAVE :: XPHI_LIM ! Threshold value for Phi3 and Psi3
REAL,SAVE :: XSBL_O_BL ! SBL height / BL height ratio
REAL,SAVE :: XFTOP_O_FSURF! Fraction of surface (heat or momentum) flux used to define top of BL
LOGICAL,SAVE :: LHARAT ! SWITCH HARATU
+LOGICAL,SAVE :: LSTATNW ! SWITCH LSTATNW
!
END MODULE MODD_CTURB
diff --git a/src/common/turb/mode_compute_function_thermo_mf.F90 b/src/common/turb/mode_compute_function_thermo_mf.F90
index 64cc93462a8539820657547cf78cff60ce859fb9..92d46fae9efe4fad66545ef8beaf97df70e8b130 100644
--- a/src/common/turb/mode_compute_function_thermo_mf.F90
+++ b/src/common/turb/mode_compute_function_thermo_mf.F90
@@ -1,6 +1,6 @@
!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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
! ######spl
MODULE MODE_COMPUTE_FUNCTION_THERMO_MF
@@ -14,14 +14,14 @@ CONTAINS
! #################################################################
!
!!
-!!**** *COMPUTE_FUNCTION_THERMO_MF* -
+!!**** *COMPUTE_FUNCTION_THERMO_MF* -
!!
!! PURPOSE
!! -------
!!
!!** METHOD
!! ------
-!!
+!!
!!
!! EXTERNAL
!! --------
@@ -35,7 +35,7 @@ CONTAINS
!!
!! AUTHOR
!! ------
-!!
+!!
!! JP Pinty *LA*
!!
!! MODIFICATIONS
@@ -44,13 +44,15 @@ CONTAINS
!! Externalisation of computations done in TURB and MF_TURB (Malardel and Pergaud, fev. 2007)
!! Optimization : V.Masson, 09/2010
!! S. Riette Sept 2011 : remove of unused PL?OCPEXN, use of received ice fraction
+!! Wim de Rooy June 2019: update statistical cloud scheme
!!
!! --------------------------------------------------------------------------
-!
+!
!* 0. DECLARATIONS
! ------------
!
USE MODD_CST
+USE MODD_CTURB, ONLY : LSTATNW
USE PARKIND1, ONLY : JPRB
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
!
@@ -73,12 +75,12 @@ REAL, DIMENSION(:,:), INTENT(OUT) :: PT ! temperature
REAL, DIMENSION(:,:), INTENT(OUT) :: PAMOIST,PATHETA
!
!-------------------------------------------------------------------------------
-!
+!
!* 0.2 Declarations of local variables
!
REAL :: ZEPS ! XMV / XMD
REAL, DIMENSION(SIZE(PTH,1),SIZE(PTH,2)) :: &
- ZCP, & ! Cp
+ ZCP, & ! Cp
ZE, & ! Saturation mixing ratio
ZDEDT, & ! Saturation mixing ratio derivative
ZAMOIST_W, & ! Coefficients for s = f (Thetal,Rnp)
@@ -103,11 +105,11 @@ ZCP=XCPD
IF (KRR > 0) ZCP(:,:) = ZCP(:,:) + XCPV * PR(:,:,1)
-DO JRR = 2,1+KRRL ! loop on the liquid components
+DO JRR = 2,1+KRRL ! loop on the liquid components
ZCP(:,:) = ZCP(:,:) + XCL * PR(:,:,JRR)
END DO
-DO JRR = 2+KRRL,1+KRRL+KRRI ! loop on the solid components
+DO JRR = 2+KRRL,1+KRRL+KRRI ! loop on the solid components
ZCP(:,:) = ZCP(:,:) + XCI * PR(:,:,JRR)
END DO
@@ -118,9 +120,9 @@ PT(:,:) = PTH(:,:) * PEXN(:,:)
!
!! Liquid water
!
-IF ( KRRL >= 1 ) THEN
+IF ( KRRL >= 1 ) THEN
!
-!* Lv/Cph
+!* Lv/Cph
!
ZLVOCP(:,:) = (XLVTT + (XCPV-XCL) * (PT(:,:)-XTT) ) / ZCP(:,:)
!
@@ -137,31 +139,32 @@ IF ( KRRL >= 1 ) THEN
ZDEDT(:,:) = ( XBETAW / PT(:,:) - XGAMW ) / PT(:,:) &
* ZE(:,:) * ( 1. + ZE(:,:) / ZEPS )
!
-!* Compute Amoist
-!
- ZAMOIST_W(:,:)= 0.5 / ( 1.0 + ZDEDT(:,:) * ZLVOCP(:,:) )
-!
-!* Compute Atheta
+!* Compute Amoist and Atheta
!
- ZATHETA_W(:,:)= ZAMOIST_W(:,:) * PEXN(:,:) * &
- ( ( ZE(:,:) - PR(:,:,1) ) * ZLVOCP(:,:) / &
- ( 1. + ZDEDT(:,:) * ZLVOCP(:,:) ) * &
- ( &
- ZE(:,:) * (1. + ZE(:,:)/ZEPS) &
- * ( -2.*XBETAW/PT(:,:) + XGAMW ) / PT(:,:)**2 &
- +ZDEDT(:,:) * (1. + 2. * ZE(:,:)/ZEPS) &
- * ( XBETAW/PT(:,:) - XGAMW ) / PT(:,:) &
- ) &
- - ZDEDT(:,:) &
- )
-
+ IF (LSTATNW) THEN
+ ZAMOIST_W(:,:)= 1.0 / ( 1.0 + ZDEDT(:,:) * ZLVOCP(:,:) )
+ ZATHETA_W(:,:)= ZAMOIST_W(:,:) * PEXN(:,:) * ZDEDT(:,:)
+ ELSE
+ ZAMOIST_W(:,:)= 0.5 / ( 1.0 + ZDEDT(:,:) * ZLVOCP(:,:) )
+ ZATHETA_W(:,:)= ZAMOIST_W(:,:) * PEXN(:,:) * &
+ ( ( ZE(:,:) - PR(:,:,1) ) * ZLVOCP(:,:) / &
+ ( 1. + ZDEDT(:,:) * ZLVOCP(:,:) ) * &
+ ( &
+ ZE(:,:) * (1. + ZE(:,:)/ZEPS) &
+ * ( -2.*XBETAW/PT(:,:) + XGAMW ) / PT(:,:)**2 &
+ +ZDEDT(:,:) * (1. + 2. * ZE(:,:)/ZEPS) &
+ * ( XBETAW/PT(:,:) - XGAMW ) / PT(:,:) &
+ ) &
+ - ZDEDT(:,:) &
+ )
+ ENDIF
!
!! Solid water
!
- IF ( KRRI >= 1 ) THEN
+ IF ( KRRI >= 1 ) THEN
!
-!* Ls/Cph
+!* Ls/Cph
!
ZLSOCP(:,:) = (XLSTT + (XCPV-XCI) * (PT(:,:)-XTT) ) / ZCP(:,:)
!
@@ -178,23 +181,25 @@ IF ( KRRL >= 1 ) THEN
ZDEDT(:,:) = ( XBETAI / PT(:,:) - XGAMI ) / PT(:,:) &
* ZE(:,:) * ( 1. + ZE(:,:) / ZEPS )
!
-!* Compute Amoist
-!
- ZAMOIST_I(:,:)= 0.5 / ( 1.0 + ZDEDT(:,:) * ZLSOCP(:,:) )
-!
-!* Compute Atheta
-!
- ZATHETA_I(:,:)= ZAMOIST_I(:,:) * PEXN(:,:) * &
- ( ( ZE(:,:) - PR(:,:,1) ) * ZLSOCP(:,:) / &
- ( 1. + ZDEDT(:,:) * ZLSOCP(:,:) ) * &
- ( &
- ZE(:,:) * (1. + ZE(:,:)/ZEPS) &
- * ( -2.*XBETAI/PT(:,:) + XGAMI ) / PT(:,:)**2 &
- +ZDEDT(:,:) * (1. + 2. * ZE(:,:)/ZEPS) &
- * ( XBETAI/PT(:,:) - XGAMI ) / PT(:,:) &
- ) &
- - ZDEDT(:,:) &
- )
+!* Compute Amoist and Atheta
+!
+ IF (LSTATNW) THEN
+ ZAMOIST_I(:,:)= 1.0 / ( 1.0 + ZDEDT(:,:) * ZLSOCP(:,:) )
+ ZATHETA_I(:,:)= ZAMOIST_I(:,:) * PEXN(:,:) * ZDEDT(:,:)
+ ELSE
+ ZAMOIST_I(:,:)= 0.5 / ( 1.0 + ZDEDT(:,:) * ZLSOCP(:,:) )
+ ZATHETA_I(:,:)= ZAMOIST_I(:,:) * PEXN(:,:) * &
+ ( ( ZE(:,:) - PR(:,:,1) ) * ZLSOCP(:,:) / &
+ ( 1. + ZDEDT(:,:) * ZLSOCP(:,:) ) * &
+ ( &
+ ZE(:,:) * (1. + ZE(:,:)/ZEPS) &
+ * ( -2.*XBETAI/PT(:,:) + XGAMI ) / PT(:,:)**2 &
+ +ZDEDT(:,:) * (1. + 2. * ZE(:,:)/ZEPS) &
+ * ( XBETAI/PT(:,:) - XGAMI ) / PT(:,:) &
+ ) &
+ - ZDEDT(:,:) &
+ )
+ ENDIF
ELSE
ZAMOIST_I(:,:)=0.
diff --git a/src/common/turb/mode_compute_mf_cloud_stat.F90 b/src/common/turb/mode_compute_mf_cloud_stat.F90
index 12fcce462fef6b0eb2c4b73fec8b6058a78ada08..7f670d5b6e4650225975d52650c648ffcff41287 100644
--- a/src/common/turb/mode_compute_mf_cloud_stat.F90
+++ b/src/common/turb/mode_compute_mf_cloud_stat.F90
@@ -48,10 +48,13 @@ CONTAINS
!! -------------
!! Original 25 Aug 2011
!! S. Riette Jan 2012: support for both order of vertical levels
+!! Wim de Rooy June 2019: update statistical cloud scheme (now including
+!! covariance term for MF contribution)
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
+USE MODD_CTURB, ONLY : LSTATNW, XCTV
USE MODD_PARAM_MFSHALL_n, ONLY : XTAUSIGMF
USE MODD_PARAMETERS, ONLY : JPHEXT, JPVEXT
!
@@ -113,8 +116,13 @@ IF (KRRL > 0) THEN
!
!
- ZFLXZ(:,:) = -2 * XTAUSIGMF * PEMF(:,:)*(PTHL_UP(:,:)-MZM_MF(PTHLM(:,:), KKA, KKU, KKL)) * &
- GZ_M_W_MF(PTHLM(:,:),PDZZ(:,:), KKA, KKU, KKL)
+ IF (LSTATNW) THEN
+ ZFLXZ(:,:) = -2 * XCTV * XTAUSIGMF * PEMF(:,:)*(PTHL_UP(:,:)-MZM_MF(PTHLM(:,:), KKA, KKU, KKL)) * &
+ GZ_M_W_MF(PTHLM(:,:),PDZZ(:,:), KKA, KKU, KKL)
+ ELSE
+ ZFLXZ(:,:) = -2 * XTAUSIGMF * PEMF(:,:)*(PTHL_UP(:,:)-MZM_MF(PTHLM(:,:), KKA, KKU, KKL)) * &
+ GZ_M_W_MF(PTHLM(:,:),PDZZ(:,:), KKA, KKU, KKL)
+ ENDIF
!
! Avoid negative values
ZFLXZ(:,:) = MAX(0.,ZFLXZ(:,:))
@@ -129,14 +137,30 @@ IF (KRRL > 0) THEN
!
!
!
- ZFLXZ(:,:) = -2 * XTAUSIGMF * PEMF(:,:)*(PRT_UP(:,:)-MZM_MF(PRTM(:,:), KKA, KKU, KKL)) * &
- GZ_M_W_MF(PRTM(:,:),PDZZ(:,:), KKA, KKU, KKL)
+ IF (LSTATNW) THEN
+ ZFLXZ(:,:) = -2 * XCTV * XTAUSIGMF * PEMF(:,:)*(PRT_UP(:,:)-MZM_MF(PRTM(:,:), KKA, KKU, KKL)) * &
+ GZ_M_W_MF(PRTM(:,:),PDZZ(:,:), KKA, KKU, KKL)
+ ELSE
+ ZFLXZ(:,:) = -2 * XTAUSIGMF * PEMF(:,:)*(PRT_UP(:,:)-MZM_MF(PRTM(:,:), KKA, KKU, KKL)) * &
+ GZ_M_W_MF(PRTM(:,:),PDZZ(:,:), KKA, KKU, KKL)
+ ENDIF
!
! Avoid negative values
ZFLXZ(:,:) = MAX(0.,ZFLXZ(:,:))
!
PSIGMF(:,:) = PSIGMF(:,:) + ZAMOIST(:,:) **2 * MZF_MF(ZFLXZ(:,:), KKA, KKU, KKL)
+ IF (LSTATNW) THEN
+ !wc Now including convection covariance contribution in case of LSTATNW=TRUE
+ !
+ ! 1.2.2 contribution from <Rnp Thl>
+ ZFLXZ(:,:) = - XCTV * XTAUSIGMF * (PEMF(:,:)*(PRT_UP(:,:)-MZM_MF(PRTM(:,:), KKA, KKU, KKL)) * &
+ GZ_M_W_MF(PTHLM(:,:),PDZZ(:,:), KKA, KKU, KKL) + &
+ PEMF(:,:)*(PTHL_UP(:,:)-MZM_MF(PTHLM(:,:), KKA, KKU, KKL)) * &
+ GZ_M_W_MF(PRTM(:,:),PDZZ(:,:), KKA, KKU, KKL))
+
+ PSIGMF(:,:) = PSIGMF(:,:) - MIN(0.,2.*ZAMOIST(:,:)*ZATHETA(:,:)*MZF_MF(ZFLXZ(:,:), KKA, KKU, KKL))
+ ENDIF
!
! 1.3 Vertical part of Sigma_s
!
diff --git a/src/common/turb/mode_turb_ver_dyn_flux.F90 b/src/common/turb/mode_turb_ver_dyn_flux.F90
index 12915e63529d9791f97de93bbb1fe0f9cad3846b..e4c1dff95ea98c9bd16bd7e63145bed84c94763e 100644
--- a/src/common/turb/mode_turb_ver_dyn_flux.F90
+++ b/src/common/turb/mode_turb_ver_dyn_flux.F90
@@ -27,13 +27,13 @@ SUBROUTINE TURB_VER_DYN_FLUX(KKA,KKU,KKL, &
!! PURPOSE
!! -------
! The purpose of this routine is to compute the vertical turbulent
-! fluxes of the evolutive variables and give back the source
+! fluxes of the evolutive variables and give back the source
! terms to the main program. In the case of large horizontal meshes,
! the divergence of these vertical turbulent fluxes represent the whole
! effect of the turbulence but when the three-dimensionnal version of
! the turbulence scheme is activated (CTURBDIM="3DIM"), these divergences
-! are completed in the next routine TURB_HOR.
-! An arbitrary degree of implicitness has been implemented for the
+! are completed in the next routine TURB_HOR.
+! An arbitrary degree of implicitness has been implemented for the
! temporal treatment of these diffusion terms.
! The vertical boundary conditions are as follows:
! * at the bottom, the surface fluxes are prescribed at the same
@@ -41,8 +41,8 @@ SUBROUTINE TURB_VER_DYN_FLUX(KKA,KKU,KKL, &
! * at the top, the turbulent fluxes are set to 0.
! It should be noted that the condensation has been implicitely included
! in this turbulence scheme by using conservative variables and computing
-! the subgrid variance of a statistical variable s indicating the presence
-! or not of condensation in a given mesh.
+! the subgrid variance of a statistical variable s indicating the presence
+! or not of condensation in a given mesh.
!
!!** METHOD
!! ------
@@ -51,27 +51,27 @@ SUBROUTINE TURB_VER_DYN_FLUX(KKA,KKU,KKL, &
!! implicit scheme (a Crank-Nicholson type with coefficients different
!! than 0.5), which allows to vary the degree of implicitness of the
!! formulation.
-!! The different prognostic variables are treated one by one.
-!! The contributions of each turbulent fluxes are cumulated into the
-!! tendency PRvarS, and into the dynamic and thermal production of
+!! The different prognostic variables are treated one by one.
+!! The contributions of each turbulent fluxes are cumulated into the
+!! tendency PRvarS, and into the dynamic and thermal production of
!! TKE if necessary.
-!!
+!!
!! In section 2 and 3, the thermodynamical fields are considered.
!! Only the turbulent fluxes of the conservative variables
-!! (Thetal and Rnp stored in PRx(:,:,:,1)) are computed.
-!! Note that the turbulent fluxes at the vertical
+!! (Thetal and Rnp stored in PRx(:,:,:,1)) are computed.
+!! Note that the turbulent fluxes at the vertical
!! boundaries are given either by the soil scheme for the surface one
-!! ( at the same instant as the others fluxes) and equal to 0 at the
-!! top of the model. The thermal production is computed by vertically
+!! ( at the same instant as the others fluxes) and equal to 0 at the
+!! top of the model. The thermal production is computed by vertically
!! averaging the turbulent flux and multiply this flux at the mass point by
!! a function ETHETA or EMOIST, which preform the transformation from the
-!! conservative variables to the virtual potential temperature.
-!!
+!! conservative variables to the virtual potential temperature.
+!!
!! In section 4, the variance of the statistical variable
-!! s indicating presence or not of condensation, is determined in function
+!! s indicating presence or not of condensation, is determined in function
!! of the turbulent moments of the conservative variables and its
-!! squarred root is stored in PSIGS. This information will be completed in
-!! the horizontal turbulence if the turbulence dimensionality is not
+!! squarred root is stored in PSIGS. This information will be completed in
+!! the horizontal turbulence if the turbulence dimensionality is not
!! equal to "1DIM".
!!
!! In section 5, the x component of the stress tensor is computed.
@@ -82,50 +82,50 @@ SUBROUTINE TURB_VER_DYN_FLUX(KKA,KKU,KKL, &
!! j" is also parallel to the surface and in the normal direction of
!! the maximum slope
!! k" is the normal to the surface
-!! In order to prevent numerical instability, the implicit scheme has
-!! been extended to the surface flux regarding to its dependence in
-!! function of U. The dependence in function of the other components
+!! In order to prevent numerical instability, the implicit scheme has
+!! been extended to the surface flux regarding to its dependence in
+!! function of U. The dependence in function of the other components
!! introduced by the different rotations is only explicit.
-!! The turbulent fluxes are used to compute the dynamic production of
+!! The turbulent fluxes are used to compute the dynamic production of
!! TKE. For the last TKE level ( located at PDZZ(:,:,IKB)/2 from the
-!! ground), an harmonic extrapolation from the dynamic production at
+!! ground), an harmonic extrapolation from the dynamic production at
!! PDZZ(:,:,IKB) is used to avoid an evaluation of the gradient of U
!! in the surface layer.
!!
!! In section 6, the same steps are repeated but for the y direction
-!! and in section 7, a diagnostic computation of the W variance is
+!! and in section 7, a diagnostic computation of the W variance is
!! performed.
!!
-!! In section 8, the turbulent fluxes for the scalar variables are
+!! In section 8, the turbulent fluxes for the scalar variables are
!! computed by the same way as the conservative thermodynamical variables
!!
-!!
+!!
!! EXTERNAL
!! --------
-!! GX_U_M, GY_V_M, GZ_W_M : cartesian gradient operators
+!! GX_U_M, GY_V_M, GZ_W_M : cartesian gradient operators
!! GX_U_UW,GY_V_VW (X,Y,Z) represent the direction of the gradient
-!! _(M,U,...)_ represent the localization of the
+!! _(M,U,...)_ represent the localization of the
!! field to be derivated
-!! _(M,UW,...) represent the localization of the
+!! _(M,UW,...) represent the localization of the
!! field derivated
-!!
+!!
!!
!! MXM,MXF,MYM,MYF,MZM,MZF
-!! : Shuman functions (mean operators)
+!! : Shuman functions (mean operators)
!! DXF,DYF,DZF,DZM
-!! : Shuman functions (difference operators)
-!!
+!! : Shuman functions (difference operators)
+!!
!! SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
!! of a variable located at a wind point
!!
-!! FUNCTIONs ETHETA and EMOIST :
+!! FUNCTIONs ETHETA and EMOIST :
!! allows to compute:
!! - the coefficients for the turbulent correlation between
-!! any variable and the virtual potential temperature, of its
-!! correlations with the conservative potential temperature and
+!! any variable and the virtual potential temperature, of its
+!! correlations with the conservative potential temperature and
!! the humidity conservative variable:
!! ------- ------- -------
-!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
+!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
!!
!!
!! IMPLICIT ARGUMENTS
@@ -159,34 +159,34 @@ SUBROUTINE TURB_VER_DYN_FLUX(KKA,KKU,KKL, &
!! MODIFICATIONS
!! -------------
!! Original August 19, 1994
-!! Modifications: February 14, 1995 (J.Cuxart and J.Stein)
+!! Modifications: February 14, 1995 (J.Cuxart and J.Stein)
!! Doctorization and Optimization
-!! Modifications: March 21, 1995 (J.M. Carriere)
+!! Modifications: March 21, 1995 (J.M. Carriere)
!! Introduction of cloud water
-!! Modifications: June 14, 1995 (J.Cuxart and J. Stein)
+!! Modifications: June 14, 1995 (J.Cuxart and J. Stein)
!! Phi3 and Psi3 at w-point + bug in the all
-!! or nothing condens.
-!! Modifications: Sept 15, 1995 (J.Cuxart and J. Stein)
+!! or nothing condens.
+!! Modifications: Sept 15, 1995 (J.Cuxart and J. Stein)
!! Change the DP computation at the ground
-!! Modifications: October 10, 1995 (J.Cuxart and J. Stein)
+!! Modifications: October 10, 1995 (J.Cuxart and J. Stein)
!! Psi for scal var and LES tools
!! Modifications: November 10, 1995 (J. Stein)
-!! change the surface relations
+!! change the surface relations
!! Modifications: February 20, 1995 (J. Stein) optimization
-!! Modifications: May 21, 1996 (J. Stein)
-!! bug in the vertical flux of the V wind
+!! Modifications: May 21, 1996 (J. Stein)
+!! bug in the vertical flux of the V wind
!! component for explicit computation
-!! Modifications: May 21, 1996 (N. wood)
+!! Modifications: May 21, 1996 (N. wood)
!! modify the computation of the vertical
!! part or the surface tangential flux
!! Modifications: May 21, 1996 (P. Jabouille)
!! same modification in the Y direction
-!!
+!!
!! Modifications: Sept 17, 1996 (J. Stein) change the moist case by using
!! Pi instead of Piref + use Atheta and Amoist
!!
-!! Modifications: Nov 24, 1997 (V. Masson) removes the DO loops
-!! Modifications: Mar 31, 1998 (V. Masson) splits the routine TURB_VER_DYN_FLUX
+!! Modifications: Nov 24, 1997 (V. Masson) removes the DO loops
+!! Modifications: Mar 31, 1998 (V. Masson) splits the routine TURB_VER_DYN_FLUX
!! Modifications: Oct 18, 2000 (J. Stein) Bug in some computations for IKB level
!! Modifications: Oct 18, 2000 (V. Masson) LES computations + LFLAT switch
!! Nov 06, 2002 (V. Masson) LES budgets
@@ -194,12 +194,14 @@ SUBROUTINE TURB_VER_DYN_FLUX(KKA,KKU,KKL, &
!! change of YCOMMENT
!! 2012-02 Y. Seity, add possibility to run with reversed vertical levels
!! Modifications July 2015 (Wim de Rooy) LHARATU switch
-!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
!! Q. Rodier 17/01/2019 : cleaning : remove cyclic conditions on DP and ZA
+!! Modification June 2019 (Wim de Rooy) 50*MF term can be removed with
+!! inclusion of energy cascade
!! JL Redelsperger 03/2021 : Add Ocean & O-A Autocoupling LES Cases
!!--------------------------------------------------------------------------
-!
+!
!* 0. DECLARATIONS
! ------------
!
@@ -238,7 +240,7 @@ IMPLICIT NONE
!
!
!
-INTEGER, INTENT(IN) :: KKA !near ground array index
+INTEGER, INTENT(IN) :: KKA !near ground array index
INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
@@ -251,14 +253,14 @@ REAL, INTENT(IN) :: PIMPL, PEXPL ! Coef. for temporal disc.
REAL, INTENT(IN) :: PTSTEP ! Double Time Step
TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
! Metric coefficients
REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSZW ! Director Cosinus of the
! normal to the ground surface
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! altitude of flux points
-REAL, DIMENSION(:,:), INTENT(IN) :: PCOSSLOPE ! cosinus of the angle
+REAL, DIMENSION(:,:), INTENT(IN) :: PCOSSLOPE ! cosinus of the angle
! between i and the slope vector
-REAL, DIMENSION(:,:), INTENT(IN) :: PSINSLOPE ! sinus of the angle
+REAL, DIMENSION(:,:), INTENT(IN) :: PSINSLOPE ! sinus of the angle
! between i and the slope vector
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry density * grid volum
@@ -266,8 +268,8 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: MFMOIST ! moist mass flux dual sch
!
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
+REAL, DIMENSION(:,:), INTENT(IN) :: PTAU11M ! <uu> in the axes linked
+ ! to the maximum slope direction and the surface normal and the binormal
! at time t - dt
REAL, DIMENSION(:,:), INTENT(IN) :: PTAU12M ! <uv> in the same axes
REAL, DIMENSION(:,:), INTENT(IN) :: PTAU33M ! <ww> in the same axes
@@ -276,13 +278,13 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PUM,PVM,PWM, PTHLM
! Wind at t-Delta t
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM
-REAL, DIMENSION(:,:), INTENT(IN) :: PUSLOPEM ! wind component along the
+REAL, DIMENSION(:,:), INTENT(IN) :: PUSLOPEM ! wind component along the
! maximum slope direction
-REAL, DIMENSION(:,:), INTENT(IN) :: PVSLOPEM ! wind component along the
+REAL, DIMENSION(:,:), INTENT(IN) :: PVSLOPEM ! wind component along the
! direction normal to the maximum slope one
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at time t
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PWU ! momentum flux u'w'
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PWV ! momentum flux v'w'
!
@@ -300,7 +302,7 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PTP ! Thermal TKE production t
!
REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2)) :: ZDIRSINZW ! sinus of the angle
! between the normal and the vertical at the surface
-REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),1):: ZCOEFS ! coeff. for the
+REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),1):: ZCOEFS ! coeff. for the
! implicit scheme for the wind at the surface
REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),SIZE(PUM,3)) :: &
ZA, & ! under diagonal elements of the tri-diagonal matrix involved
@@ -308,7 +310,7 @@ REAL, DIMENSION(SIZE(PUM,1),SIZE(PUM,2),SIZE(PUM,3)) :: &
! J in Section 5)
ZRES, & ! guess of the treated variable at t+ deltat when the turbu-
! lence is the only source of evolution added to the ones
- ! considered in ZSOURCE
+ ! considered in ZSOURCE
ZFLXZ, & ! vertical flux of the treated variable
ZSOURCE, & ! source of evolution for the treated variable
ZKEFF ! effectif diffusion coeff = LT * SQRT( TKE )
@@ -322,7 +324,7 @@ REAL, DIMENSION(SIZE(PDZZ,1),SIZE(PDZZ,2),1) :: ZCOEFFLXU, &
ZCOEFFLXV, ZUSLOPEM, ZVSLOPEM
! coefficients for the surface flux
! evaluation and copy of PUSLOPEM and
- ! PVSLOPEM in local 3D arrays
+ ! PVSLOPEM in local 3D arrays
INTEGER :: IIU,IJU ! size of array in x,y,z directions
!
REAL :: ZTIME1, ZTIME2, ZCMFS
@@ -344,8 +346,8 @@ IKB=KKA+JPVEXT_TURB*KKL
IKE=KKU-JPVEXT_TURB*KKL
IKB=KKA+JPVEXT_TURB*KKL
IKE=KKU-JPVEXT_TURB*KKL
-IKT=SIZE(PUM,3)
-IKTB=1+JPVEXT_TURB
+IKT=SIZE(PUM,3)
+IKTB=1+JPVEXT_TURB
IKTE=IKT-JPVEXT_TURB
!
ZSOURCE = 0.
@@ -354,14 +356,16 @@ ZCMFS = XCMFS
IF (LHARAT) ZCMFS=1.
!
ZDIRSINZW(:,:) = SQRT(1.-PDIRCOSZW(:,:)**2)
-! compute the coefficients for the uncentred gradient computation near the
+! compute the coefficients for the uncentred gradient computation near the
! ground
!
! With LHARATU length scale and TKE are at half levels so remove MZM
!
IF (LHARAT) THEN
- ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
-ELSE
+ !wc with energy cascade 50*MF no longer necessary
+ !ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
+ ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:))
+ELSE
ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
ENDIF
@@ -372,21 +376,21 @@ ZVSLOPEM(:,:,1)=PVSLOPEM(:,:)
!----------------------------------------------------------------------------
!
!
-!* 5. SOURCES OF U,W WIND COMPONENTS AND PARTIAL DYNAMIC PRODUCTION
+!* 5. SOURCES OF U,W WIND COMPONENTS AND PARTIAL DYNAMIC PRODUCTION
! -------------------------------------------------------------
!
!* 5.1 Source of U wind component
!
-! Preparation of the arguments for TRIDIAG_WIND
+! Preparation of the arguments for TRIDIAG_WIND
!
ZA(:,:,:) = -PTSTEP * ZCMFS * &
MXM( ZKEFF ) * MXM(MZM(PRHODJ, KKA, KKU, KKL)) / &
MXM( PDZZ )**2
!
!
-! Compute the source of U wind component
+! Compute the source of U wind component
!
-! compute the coefficient between the vertical flux and the 2 components of the
+! compute the coefficient between the vertical flux and the 2 components of the
! wind following the slope
ZCOEFFLXU(:,:,1) = PCDUEFF(:,:) * (PDIRCOSZW(:,:)**2 - ZDIRSINZW(:,:)**2) &
* PCOSSLOPE(:,:)
@@ -403,9 +407,9 @@ ZCOEFS(:,:,1:1)=MXM(ZCOEFS(:,:,1:1) / PDZZ(:,:,IKB:IKB) )
! ZSOURCE= FLUX /DZ
IF (OOCEAN) THEN ! OCEAN MODEL ONLY
! Sfx flux assumed to be in SI & at vorticity point
- IF (LCOUPLES) THEN
+ IF (LCOUPLES) THEN
ZSOURCE(:,:,IKE:IKE) = XSSUFL_C(:,:,1:1)/PDZZ(:,:,IKE:IKE) &
- *0.5 * ( 1. + MXM(PRHODJ(:,:,KKU:KKU)) / MXM(PRHODJ(:,:,IKE:IKE)))
+ *0.5 * ( 1. + MXM(PRHODJ(:,:,KKU:KKU)) / MXM(PRHODJ(:,:,IKE:IKE)))
ELSE
ZSOURCE(:,:,IKE) = XSSUFL(:,:)
ZSOURCE(:,:,IKE:IKE) = ZSOURCE (:,:,IKE:IKE) /PDZZ(:,:,IKE:IKE) &
@@ -416,15 +420,15 @@ IF (OOCEAN) THEN ! OCEAN MODEL ONLY
ZSOURCE(:,:,IKTB+1:IKTE-1) = 0
!
ELSE !ATMOS MODEL ONLY
- IF (LCOUPLES) THEN
+ IF (LCOUPLES) THEN
ZSOURCE(:,:,IKB:IKB) = XSSUFL_C(:,:,1:1)/PDZZ(:,:,IKB:IKB) &
* 0.5 * ( 1. + MXM(PRHODJ(:,:,KKA:KKA)) / MXM(PRHODJ(:,:,IKB:IKB)) )
- ELSE
+ ELSE
! compute the explicit tangential flux at the W point
ZSOURCE(:,:,IKB) = &
PTAU11M(:,:) * PCOSSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
-PTAU12M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) &
- -PTAU33M(:,:) * PCOSSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
+ -PTAU33M(:,:) * PCOSSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
!
! add the vertical part or the surface flux at the U,W vorticity point
!
@@ -436,17 +440,17 @@ ELSE !ATMOS MODEL ONLY
*ZVSLOPEM(:,:,1:1) ) &
- ZCOEFS(:,:,1:1) * PUM(:,:,IKB:IKB) * PIMPL &
) * 0.5 * ( 1. + MXM(PRHODJ(:,:,KKA:KKA)) / MXM(PRHODJ(:,:,IKB:IKB)) )
- ENDIF
+ ENDIF
!
ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
ZSOURCE(:,:,IKE) = 0.
ENDIF !end ocean or atmosphere cases
!
-! Obtention of the split U at t+ deltat
+! Obtention of the split U at t+ deltat
!
CALL TRIDIAG_WIND(KKA,KKU,KKL,PUM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL, &
MXM(PRHODJ),ZSOURCE,ZRES)
-!
+!
! Compute the equivalent tendency for the U wind component
!
PRUS(:,:,:)=PRUS(:,:,:)+MXM(PRHODJ(:,:,:))*(ZRES(:,:,:)-PUM(:,:,:))/PTSTEP
@@ -459,19 +463,19 @@ PRUS(:,:,:)=PRUS(:,:,:)+MXM(PRHODJ(:,:,:))*(ZRES(:,:,:)-PUM(:,:,:))/PTSTEP
ZFLXZ(:,:,:) = -ZCMFS * MXM(ZKEFF) * &
DZM(PIMPL*ZRES + PEXPL*PUM, KKA, KKU, KKL) / MXM(PDZZ)
!
-! surface flux
+! surface flux
ZFLXZ(:,:,IKB:IKB) = MXM(PDZZ(:,:,IKB:IKB)) * &
( ZSOURCE(:,:,IKB:IKB) &
- +ZCOEFS(:,:,1:1) * ZRES(:,:,IKB:IKB) * PIMPL &
+ +ZCOEFS(:,:,1:1) * ZRES(:,:,IKB:IKB) * PIMPL &
) / 0.5 / ( 1. + MXM(PRHODJ(:,:,KKA:KKA)) / MXM(PRHODJ(:,:,IKB:IKB)) )
!
-ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
IF (OOCEAN) THEN !ocean model at phys sfc (ocean domain top)
ZFLXZ(:,:,IKE:IKE) = MXM(PDZZ(:,:,IKE:IKE)) * &
ZSOURCE(:,:,IKE:IKE) &
/ 0.5 / ( 1. + MXM(PRHODJ(:,:,KKU:KKU)) / MXM(PRHODJ(:,:,IKE:IKE)) )
- ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE)
+ ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE)
END IF
!
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
@@ -509,14 +513,14 @@ IF (OOCEAN) THEN
PDP(:,:,IKE:IKE) = - MXF ( &
ZFLXZ(:,:,IKE-KKL:IKE-KKL) * (PUM(:,:,IKE:IKE)-PUM(:,:,IKE-KKL:IKE-KKL)) &
/ MXM(PDZZ(:,:,IKE-KKL:IKE-KKL)) &
- )
+ )
END IF
!
! Storage in the LES configuration
-!
+!
IF (LLES_CALL) THEN
CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID(MZF(MXF(ZFLXZ), KKA, KKU, KKL), X_LES_SUBGRID_WU )
+ CALL LES_MEAN_SUBGRID(MZF(MXF(ZFLXZ), KKA, KKU, KKL), X_LES_SUBGRID_WU )
CALL LES_MEAN_SUBGRID(MZF(MXF(GZ_U_UW(PUM,PDZZ, KKA, KKU, KKL) &
& *ZFLXZ), KKA, KKU, KKL), X_LES_RES_ddxa_U_SBG_UaU )
CALL LES_MEAN_SUBGRID( ZCMFS * ZKEFF, X_LES_SUBGRID_Km )
@@ -530,11 +534,11 @@ END IF
IF(HTURBDIM=='3DIM') THEN
! Compute the source for the W wind component
! used to compute the W source at the ground
- ZFLXZ(:,:,KKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation
+ ZFLXZ(:,:,KKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation
IF (OOCEAN) THEN
- ZFLXZ(:,:,KKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-KKL) ! extrapolation
- END IF
-
+ ZFLXZ(:,:,KKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-KKL) ! extrapolation
+ END IF
+
!
IF (.NOT. LFLAT) THEN
PRWS(:,:,:)= PRWS &
@@ -546,7 +550,7 @@ IF(HTURBDIM=='3DIM') THEN
PRWS(:,:,:)= PRWS -DXF(MZM(MXM(PRHODJ) /PDXX, KKA, KKU, KKL) * ZFLXZ )
END IF
!
- ! Complete the Dynamical production with the W wind component
+ ! Complete the Dynamical production with the W wind component
!
ZA(:,:,:)=-MZF(MXF(ZFLXZ * GX_W_UW(PWM,PDXX,PDZZ,PDZX, KKA, KKU, KKL)), KKA, KKU, KKL)
!
@@ -582,7 +586,7 @@ END IF
PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
!
! Storage in the LES configuration
- !
+ !
IF (LLES_CALL) THEN
CALL SECOND_MNH(ZTIME1)
CALL LES_MEAN_SUBGRID(MZF(MXF(GX_W_UW(PWM,PDXX,&
@@ -605,7 +609,7 @@ END IF
!----------------------------------------------------------------------------
!
!
-!* 6. SOURCES OF V,W WIND COMPONENTS AND COMPLETE 1D DYNAMIC PRODUCTION
+!* 6. SOURCES OF V,W WIND COMPONENTS AND COMPLETE 1D DYNAMIC PRODUCTION
! -----------------------------------------------------------------
!
!* 6.1 Source of V wind component
@@ -619,7 +623,7 @@ ZA(:,:,:) = - PTSTEP * ZCMFS * &
!
!
! Compute the source of V wind component
-! compute the coefficient between the vertical flux and the 2 components of the
+! compute the coefficient between the vertical flux and the 2 components of the
! wind following the slope
ZCOEFFLXU(:,:,1) = PCDUEFF(:,:) * (PDIRCOSZW(:,:)**2 - ZDIRSINZW(:,:)**2) &
* PSINSLOPE(:,:)
@@ -635,8 +639,8 @@ ZCOEFS(:,:,1:1)=MYM(ZCOEFS(:,:,1:1) / PDZZ(:,:,IKB:IKB) )
IF (OOCEAN) THEN ! Ocean case
IF (LCOUPLES) THEN
ZSOURCE(:,:,IKE:IKE) = XSSVFL_C(:,:,1:1)/PDZZ(:,:,IKE:IKE) &
- *0.5 * ( 1. + MYM(PRHODJ(:,:,KKU:KKU)) / MYM(PRHODJ(:,:,IKE:IKE)) )
- ELSE
+ *0.5 * ( 1. + MYM(PRHODJ(:,:,KKU:KKU)) / MYM(PRHODJ(:,:,IKE:IKE)) )
+ ELSE
ZSOURCE(:,:,IKE) = XSSVFL(:,:)
ZSOURCE(:,:,IKE:IKE) = ZSOURCE(:,:,IKE:IKE)/PDZZ(:,:,IKE:IKE) &
*0.5 * ( 1. + MYM(PRHODJ(:,:,KKU:KKU)) / MYM(PRHODJ(:,:,IKE:IKE)) )
@@ -649,7 +653,7 @@ ELSE ! Atmos case
ZSOURCE(:,:,IKB) = &
PTAU11M(:,:) * PSINSLOPE(:,:) * PDIRCOSZW(:,:) * ZDIRSINZW(:,:) &
+PTAU12M(:,:) * PCOSSLOPE(:,:) * ZDIRSINZW(:,:) &
- -PTAU33M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
+ -PTAU33M(:,:) * PSINSLOPE(:,:) * ZDIRSINZW(:,:) * PDIRCOSZW(:,:)
!
! add the vertical part or the surface flux at the V,W vorticity point
ZSOURCE(:,:,IKB:IKB) = &
@@ -670,8 +674,8 @@ ELSE ! Atmos case
ZSOURCE(:,:,IKE) = 0.
ENDIF ! End of Ocean or Atmospher Cases
ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
-!
-! Obtention of the split V at t+ deltat
+!
+! Obtention of the split V at t+ deltat
CALL TRIDIAG_WIND(KKA,KKU,KKL,PVM,ZA,ZCOEFS(:,:,1),PTSTEP,PEXPL,PIMPL, &
MYM(PRHODJ),ZSOURCE,ZRES)
!
@@ -689,9 +693,9 @@ ZFLXZ(:,:,:) = -ZCMFS * MYM(ZKEFF) * &
!
ZFLXZ(:,:,IKB:IKB) = MYM(PDZZ(:,:,IKB:IKB)) * &
( ZSOURCE(:,:,IKB:IKB) &
- +ZCOEFS(:,:,1:1) * ZRES(:,:,IKB:IKB) * PIMPL &
+ +ZCOEFS(:,:,1:1) * ZRES(:,:,IKB:IKB) * PIMPL &
) / 0.5 / ( 1. + MYM(PRHODJ(:,:,KKA:KKA)) / MYM(PRHODJ(:,:,IKB:IKB)) )
-!
+!
!
ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
!
@@ -699,7 +703,7 @@ IF (OOCEAN) THEN
ZFLXZ(:,:,IKE:IKE) = MYM(PDZZ(:,:,IKE:IKE)) * &
ZSOURCE(:,:,IKE:IKE) &
/ 0.5 / ( 1. + MYM(PRHODJ(:,:,KKU:KKU)) / MYM(PRHODJ(:,:,IKE:IKE)) )
- ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE)
+ ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE)
END IF
!
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
@@ -747,7 +751,7 @@ PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
!
IF (LLES_CALL) THEN
CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID(MZF(MYF(ZFLXZ), KKA, KKU, KKL), X_LES_SUBGRID_WV )
+ CALL LES_MEAN_SUBGRID(MZF(MYF(ZFLXZ), KKA, KKU, KKL), X_LES_SUBGRID_WV )
CALL LES_MEAN_SUBGRID(MZF(MYF(GZ_V_VW(PVM,PDZZ, KKA, KKU, KKL)*&
& ZFLXZ), KKA, KKU, KKL), X_LES_RES_ddxa_V_SBG_UaV )
CALL SECOND_MNH(ZTIME2)
@@ -755,16 +759,16 @@ IF (LLES_CALL) THEN
END IF
!
!
-!* 6.3 Source of W wind component
+!* 6.3 Source of W wind component
!
IF(HTURBDIM=='3DIM') THEN
! Compute the source for the W wind component
- ZFLXZ(:,:,KKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation
+ ZFLXZ(:,:,KKA) = 2 * ZFLXZ(:,:,IKB) - ZFLXZ(:,:,IKB+KKL) ! extrapolation
IF (OOCEAN) THEN
- ZFLXZ(:,:,KKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-KKL) ! extrapolation
+ ZFLXZ(:,:,KKU) = 2 * ZFLXZ(:,:,IKE) - ZFLXZ(:,:,IKE-KKL) ! extrapolation
END IF
!
- IF (.NOT. L2D) THEN
+ IF (.NOT. L2D) THEN
IF (.NOT. LFLAT) THEN
PRWS(:,:,:)= PRWS(:,:,:) &
-DYF( MZM(MYM(PRHODJ) /PDYY, KKA, KKU, KKL) * ZFLXZ ) &
@@ -775,8 +779,8 @@ IF(HTURBDIM=='3DIM') THEN
PRWS(:,:,:)= PRWS(:,:,:) -DYF(MZM(MYM(PRHODJ) /PDYY, KKA, KKU, KKL) * ZFLXZ )
END IF
END IF
- !
- ! Complete the Dynamical production with the W wind component
+ !
+ ! Complete the Dynamical production with the W wind component
IF (.NOT. L2D) THEN
ZA(:,:,:) = - MZF(MYF(ZFLXZ * GY_W_VW(PWM,PDYY,PDZZ,PDZY, KKA, KKU, KKL)), KKA, KKU, KKL)
!
@@ -806,7 +810,7 @@ IF(HTURBDIM=='3DIM') THEN
) / (0.5*(PDYY(:,:,IKE-KKL:IKE-KKL)+PDYY(:,:,IKE:IKE))) &
)
END IF
-!
+!
PDP(:,:,:)=PDP(:,:,:)+ZA(:,:,:)
!
END IF
@@ -842,7 +846,7 @@ IF ( OTURB_FLX .AND. TPFILE%LOPENED .AND. HTURBDIM == '1DIM') THEN
ZFLXZ(:,:,:)= (2./3.) * PTKEM(:,:,:) &
-ZCMFS*PLM(:,:,:)*SQRT(PTKEM(:,:,:))*GZ_W_M(PWM,PDZZ, KKA, KKU, KKL)
! to be tested &
- ! +XCMFB*(4./3.)*PLM(:,:,:)/SQRT(PTKEM(:,:,:))*PTP(:,:,:)
+ ! +XCMFB*(4./3.)*PLM(:,:,:)/SQRT(PTKEM(:,:,:))*PTP(:,:,:)
! stores the W variance
TZFIELD%CMNHNAME = 'W_VVAR'
TZFIELD%CSTDNAME = ''
diff --git a/src/common/turb/mode_turb_ver_sv_flux.F90 b/src/common/turb/mode_turb_ver_sv_flux.F90
index 271ee734e6aabd3d2cda09b7eba343113fa6e402..e7c47ec8323c67ce0b6e3352e788b450bcba1fda 100644
--- a/src/common/turb/mode_turb_ver_sv_flux.F90
+++ b/src/common/turb/mode_turb_ver_sv_flux.F90
@@ -25,13 +25,13 @@ SUBROUTINE TURB_VER_SV_FLUX(KKA,KKU,KKL, &
!! PURPOSE
!! -------
! The purpose of this routine is to compute the vertical turbulent
-! fluxes of the evolutive variables and give back the source
+! fluxes of the evolutive variables and give back the source
! terms to the main program. In the case of large horizontal meshes,
! the divergence of these vertical turbulent fluxes represent the whole
! effect of the turbulence but when the three-dimensionnal version of
! the turbulence scheme is activated (CTURBDIM="3DIM"), these divergences
-! are completed in the next routine TURB_HOR.
-! An arbitrary degree of implicitness has been implemented for the
+! are completed in the next routine TURB_HOR.
+! An arbitrary degree of implicitness has been implemented for the
! temporal treatment of these diffusion terms.
! The vertical boundary conditions are as follows:
! * at the bottom, the surface fluxes are prescribed at the same
@@ -39,8 +39,8 @@ SUBROUTINE TURB_VER_SV_FLUX(KKA,KKU,KKL, &
! * at the top, the turbulent fluxes are set to 0.
! It should be noted that the condensation has been implicitely included
! in this turbulence scheme by using conservative variables and computing
-! the subgrid variance of a statistical variable s indicating the presence
-! or not of condensation in a given mesh.
+! the subgrid variance of a statistical variable s indicating the presence
+! or not of condensation in a given mesh.
!
!!** METHOD
!! ------
@@ -49,27 +49,27 @@ SUBROUTINE TURB_VER_SV_FLUX(KKA,KKU,KKL, &
!! implicit scheme (a Crank-Nicholson type with coefficients different
!! than 0.5), which allows to vary the degree of implicitness of the
!! formulation.
-!! The different prognostic variables are treated one by one.
-!! The contributions of each turbulent fluxes are cumulated into the
-!! tendency PRvarS, and into the dynamic and thermal production of
+!! The different prognostic variables are treated one by one.
+!! The contributions of each turbulent fluxes are cumulated into the
+!! tendency PRvarS, and into the dynamic and thermal production of
!! TKE if necessary.
-!!
+!!
!! In section 2 and 3, the thermodynamical fields are considered.
!! Only the turbulent fluxes of the conservative variables
-!! (Thetal and Rnp stored in PRx(:,:,:,1)) are computed.
-!! Note that the turbulent fluxes at the vertical
+!! (Thetal and Rnp stored in PRx(:,:,:,1)) are computed.
+!! Note that the turbulent fluxes at the vertical
!! boundaries are given either by the soil scheme for the surface one
-!! ( at the same instant as the others fluxes) and equal to 0 at the
-!! top of the model. The thermal production is computed by vertically
+!! ( at the same instant as the others fluxes) and equal to 0 at the
+!! top of the model. The thermal production is computed by vertically
!! averaging the turbulent flux and multiply this flux at the mass point by
!! a function ETHETA or EMOIST, which preform the transformation from the
-!! conservative variables to the virtual potential temperature.
-!!
+!! conservative variables to the virtual potential temperature.
+!!
!! In section 4, the variance of the statistical variable
-!! s indicating presence or not of condensation, is determined in function
+!! s indicating presence or not of condensation, is determined in function
!! of the turbulent moments of the conservative variables and its
-!! squarred root is stored in PSIGS. This information will be completed in
-!! the horizontal turbulence if the turbulence dimensionality is not
+!! squarred root is stored in PSIGS. This information will be completed in
+!! the horizontal turbulence if the turbulence dimensionality is not
!! equal to "1DIM".
!!
!! In section 5, the x component of the stress tensor is computed.
@@ -80,53 +80,53 @@ SUBROUTINE TURB_VER_SV_FLUX(KKA,KKU,KKL, &
!! j" is also parallel to the surface and in the normal direction of
!! the maximum slope
!! k" is the normal to the surface
-!! In order to prevent numerical instability, the implicit scheme has
-!! been extended to the surface flux regarding to its dependence in
-!! function of U. The dependence in function of the other components
+!! In order to prevent numerical instability, the implicit scheme has
+!! been extended to the surface flux regarding to its dependence in
+!! function of U. The dependence in function of the other components
!! introduced by the different rotations is only explicit.
-!! The turbulent fluxes are used to compute the dynamic production of
+!! The turbulent fluxes are used to compute the dynamic production of
!! TKE. For the last TKE level ( located at PDZZ(:,:,IKB)/2 from the
-!! ground), an harmonic extrapolation from the dynamic production at
+!! ground), an harmonic extrapolation from the dynamic production at
!! PDZZ(:,:,IKB) is used to avoid an evaluation of the gradient of U
!! in the surface layer.
!!
!! In section 6, the same steps are repeated but for the y direction
-!! and in section 7, a diagnostic computation of the W variance is
+!! and in section 7, a diagnostic computation of the W variance is
!! performed.
!!
-!! In section 8, the turbulent fluxes for the scalar variables are
+!! In section 8, the turbulent fluxes for the scalar variables are
!! computed by the same way as the conservative thermodynamical variables
!!
-!!
+!!
!! EXTERNAL
!! --------
-!! GX_U_M, GY_V_M, GZ_W_M : cartesian gradient operators
+!! GX_U_M, GY_V_M, GZ_W_M : cartesian gradient operators
!! GX_U_UW,GY_V_VW (X,Y,Z) represent the direction of the gradient
-!! _(M,U,...)_ represent the localization of the
+!! _(M,U,...)_ represent the localization of the
!! field to be derivated
-!! _(M,UW,...) represent the localization of the
+!! _(M,UW,...) represent the localization of the
!! field derivated
-!!
+!!
!!
!! MXM,MXF,MYM,MYF,MZM,MZF
-!! : Shuman functions (mean operators)
+!! : Shuman functions (mean operators)
!! DXF,DYF,DZF,DZM
-!! : Shuman functions (difference operators)
-!!
+!! : Shuman functions (difference operators)
+!!
!! SUBROUTINE TRIDIAG : to compute the split implicit evolution
!! of a variable located at a mass point
!!
!! SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
!! of a variable located at a wind point
!!
-!! FUNCTIONs ETHETA and EMOIST :
+!! FUNCTIONs ETHETA and EMOIST :
!! allows to compute:
!! - the coefficients for the turbulent correlation between
-!! any variable and the virtual potential temperature, of its
-!! correlations with the conservative potential temperature and
+!! any variable and the virtual potential temperature, of its
+!! correlations with the conservative potential temperature and
!! the humidity conservative variable:
!! ------- ------- -------
-!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
+!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
!!
!!
!! IMPLICIT ARGUMENTS
@@ -160,34 +160,34 @@ SUBROUTINE TURB_VER_SV_FLUX(KKA,KKU,KKL, &
!! MODIFICATIONS
!! -------------
!! Original August 19, 1994
-!! Modifications: February 14, 1995 (J.Cuxart and J.Stein)
+!! Modifications: February 14, 1995 (J.Cuxart and J.Stein)
!! Doctorization and Optimization
-!! Modifications: March 21, 1995 (J.M. Carriere)
+!! Modifications: March 21, 1995 (J.M. Carriere)
!! Introduction of cloud water
-!! Modifications: June 14, 1995 (J.Cuxart and J. Stein)
+!! Modifications: June 14, 1995 (J.Cuxart and J. Stein)
!! Phi3 and Psi3 at w-point + bug in the all
-!! or nothing condens.
-!! Modifications: Sept 15, 1995 (J.Cuxart and J. Stein)
+!! or nothing condens.
+!! Modifications: Sept 15, 1995 (J.Cuxart and J. Stein)
!! Change the DP computation at the ground
-!! Modifications: October 10, 1995 (J.Cuxart and J. Stein)
+!! Modifications: October 10, 1995 (J.Cuxart and J. Stein)
!! Psi for scal var and LES tools
!! Modifications: November 10, 1995 (J. Stein)
-!! change the surface relations
+!! change the surface relations
!! Modifications: February 20, 1995 (J. Stein) optimization
-!! Modifications: May 21, 1996 (J. Stein)
-!! bug in the vertical flux of the V wind
+!! Modifications: May 21, 1996 (J. Stein)
+!! bug in the vertical flux of the V wind
!! component for explicit computation
-!! Modifications: May 21, 1996 (N. wood)
+!! Modifications: May 21, 1996 (N. wood)
!! modify the computation of the vertical
!! part or the surface tangential flux
!! Modifications: May 21, 1996 (P. Jabouille)
!! same modification in the Y direction
-!!
+!!
!! Modifications: Sept 17, 1996 (J. Stein) change the moist case by using
!! Pi instead of Piref + use Atheta and Amoist
!!
!! Modifications: Nov 24, 1997 (V. Masson) removes the DO loops
-!! Modifications: Mar 31, 1998 (V. Masson) splits the routine TURB_VER_SV_FLUX
+!! Modifications: Mar 31, 1998 (V. Masson) splits the routine TURB_VER_SV_FLUX
!! Modifications: Dec 01, 2000 (V. Masson) conservation of scalar emission
!! from surface in 1DIM case
!! when slopes are present
@@ -195,15 +195,17 @@ SUBROUTINE TURB_VER_SV_FLUX(KKA,KKU,KKL, &
!! Nov 06, 2002 (V. Masson) LES budgets
!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
!! change of YCOMMENT
-!! Feb 2012(Y. Seity) add possibility to run with reversed
+!! Feb 2012(Y. Seity) add possibility to run with reversed
!! vertical levels
!! Modifications: July 2015 (Wim de Rooy) LHARAT switch
!! Feb 2017(M. Leriche) add initialisation of ZSOURCE
!! to avoid unknwon values outside physical domain
!! and avoid negative values in sv tendencies
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Modifications: June 2019 (Wim de Rooy) with energycascade, 50MF nog
+!! longer necessary
!!--------------------------------------------------------------------------
-!
+!
!* 0. DECLARATIONS
! ------------
!
@@ -238,7 +240,7 @@ IMPLICIT NONE
!* 0.1 declarations of arguments
!
!
-INTEGER, INTENT(IN) :: KKA !near ground array index
+INTEGER, INTENT(IN) :: KKA !near ground array index
INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
@@ -258,15 +260,15 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry density * grid volum
REAL, DIMENSION(:,:,:), INTENT(IN) :: MFMOIST ! moist mf dual scheme
!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PSFSVM ! t - deltat
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSFSVM ! t - deltat
!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PSFSVP ! t + deltat
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSFSVP ! t + deltat
!
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! scalar var. at t-Delta t
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM ! vertical wind
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at time t
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PPSI_SV ! Inv.Turb.Sch.for scalars
!
REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRSVS
@@ -285,14 +287,14 @@ REAL, DIMENSION(SIZE(PSVM,1),SIZE(PSVM,2),SIZE(PSVM,3)) :: &
! J in Section 5)
ZRES, & ! guess of the treated variable at t+ deltat when the turbu-
! lence is the only source of evolution added to the ones
- ! considered in ZSOURCE
+ ! considered in ZSOURCE
ZFLXZ, & ! vertical flux of the treated variable
ZSOURCE, & ! source of evolution for the treated variable
ZKEFF ! effectif diffusion coeff = LT * SQRT( TKE )
INTEGER :: IKB,IKE ! I index values for the Beginning and End
! mass points of the domain in the 3 direct.
INTEGER :: IKT ! array size in k direction
-INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
+INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
INTEGER :: JSV ! loop counters
INTEGER :: JK ! loop
INTEGER :: ISV ! number of scalar var.
@@ -314,19 +316,21 @@ IF (LHOOK) CALL DR_HOOK('TURB_VER_SV_FLUX',0,ZHOOK_HANDLE)
IKB=KKA+JPVEXT_TURB*KKL
IKE=KKU-JPVEXT_TURB*KKL
IKT=SIZE(PSVM,3)
-IKTE =IKT-JPVEXT_TURB
-IKTB =1+JPVEXT_TURB
+IKTE =IKT-JPVEXT_TURB
+IKTB =1+JPVEXT_TURB
!
ISV=SIZE(PSVM,4)
!
IF (LHARAT) THEN
- ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50.*MFMOIST(:,:,:)
+ !wc 50MF can be omitted with energy cascade included
+ !ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50.*MFMOIST(:,:,:)
+ ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:))
ELSE
ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
ENDIF
!
IF(LBLOWSNOW) THEN
-! See Vionnet (PhD, 2012) for a complete discussion around the value of the Schmidt number for blowing snow variables
+! 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
@@ -340,7 +344,7 @@ DO JSV=1,ISV
!
IF (LNOMIXLG .AND. JSV >= NSV_LGBEG .AND. JSV<= NSV_LGEND) CYCLE
!
-! Preparation of the arguments for TRIDIAG
+! Preparation of the arguments for TRIDIAG
IF (LHARAT) THEN
ZA(:,:,:) = -PTSTEP* &
ZKEFF * MZM(PRHODJ, KKA, KKU, KKL) / &
@@ -361,7 +365,7 @@ DO JSV=1,ISV
IF (HTURBDIM=='3DIM') THEN
ZSOURCE(:,:,IKB) = (PIMPL*PSFSVP(:,:,JSV) + PEXPL*PSFSVM(:,:,JSV)) / &
PDZZ(:,:,IKB) * PDIRCOSZW(:,:) &
- * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
+ * 0.5 * (1. + PRHODJ(:,:,KKA) / PRHODJ(:,:,IKB))
ELSE
ZSOURCE(:,:,IKB) = (PIMPL*PSFSVP(:,:,JSV) + PEXPL*PSFSVM(:,:,JSV)) / &
@@ -371,7 +375,7 @@ DO JSV=1,ISV
ZSOURCE(:,:,IKTB+1:IKTE-1) = 0.
ZSOURCE(:,:,IKE) = 0.
!
-! Obtention of the split JSV scalar variable at t+ deltat
+! Obtention of the split JSV scalar variable at t+ deltat
CALL TRIDIAG(KKA,KKU,KKL,PSVM(:,:,:,JSV),ZA,PTSTEP,PEXPL,PIMPL,PRHODJ,ZSOURCE,ZRES)
!
! Compute the equivalent tendency for the JSV scalar variable
@@ -390,7 +394,7 @@ DO JSV=1,ISV
! is taken into account in the vertical part
IF (HTURBDIM=='3DIM') THEN
ZFLXZ(:,:,IKB) = (PIMPL*PSFSVP(:,:,JSV) + PEXPL*PSFSVM(:,:,JSV)) &
- * PDIRCOSZW(:,:)
+ * PDIRCOSZW(:,:)
ELSE
ZFLXZ(:,:,IKB) = (PIMPL*PSFSVP(:,:,JSV) + PEXPL*PSFSVM(:,:,JSV)) &
/ PDIRCOSZW(:,:)
@@ -437,7 +441,7 @@ DO JSV=1,ISV
XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
END IF
!
-END DO ! end of scalar loop
+END DO ! end of scalar loop
!
!----------------------------------------------------------------------------
!
diff --git a/src/common/turb/mode_turb_ver_thermo_corr.F90 b/src/common/turb/mode_turb_ver_thermo_corr.F90
index db08f0ce62947783ab09870fd2b49ff4d4245f2d..596c8c7b3af1c5b37fdd8759aa026b5aaf9cd99f 100644
--- a/src/common/turb/mode_turb_ver_thermo_corr.F90
+++ b/src/common/turb/mode_turb_ver_thermo_corr.F90
@@ -4,7 +4,7 @@
!MNH_LIC for details. version 1.
MODULE MODE_TURB_VER_THERMO_CORR
IMPLICIT NONE
-CONTAINS
+CONTAINS
SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI, &
OTURB_FLX,HTURBDIM,HTOM, &
PIMPL,PEXPL, &
@@ -29,13 +29,13 @@ SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI, &
!! PURPOSE
!! -------
! The purpose of this routine is to compute the vertical turbulent
-! fluxes of the evolutive variables and give back the source
+! fluxes of the evolutive variables and give back the source
! terms to the main program. In the case of large horizontal meshes,
! the divergence of these vertical turbulent fluxes represent the whole
! effect of the turbulence but when the three-dimensionnal version of
! the turbulence scheme is activated (CTURBDIM="3DIM"), these divergences
-! are completed in the next routine TURB_HOR.
-! An arbitrary degree of implicitness has been implemented for the
+! are completed in the next routine TURB_HOR.
+! An arbitrary degree of implicitness has been implemented for the
! temporal treatment of these diffusion terms.
! The vertical boundary conditions are as follows:
! * at the bottom, the surface fluxes are prescribed at the same
@@ -43,8 +43,8 @@ SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI, &
! * at the top, the turbulent fluxes are set to 0.
! It should be noted that the condensation has been implicitely included
! in this turbulence scheme by using conservative variables and computing
-! the subgrid variance of a statistical variable s indicating the presence
-! or not of condensation in a given mesh.
+! the subgrid variance of a statistical variable s indicating the presence
+! or not of condensation in a given mesh.
!
!!** METHOD
!! ------
@@ -53,27 +53,27 @@ SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI, &
!! implicit scheme (a Crank-Nicholson type with coefficients different
!! than 0.5), which allows to vary the degree of implicitness of the
!! formulation.
-!! The different prognostic variables are treated one by one.
-!! The contributions of each turbulent fluxes are cumulated into the
-!! tendency PRvarS, and into the dynamic and thermal production of
+!! The different prognostic variables are treated one by one.
+!! The contributions of each turbulent fluxes are cumulated into the
+!! tendency PRvarS, and into the dynamic and thermal production of
!! TKE if necessary.
-!!
+!!
!! In section 2 and 3, the thermodynamical fields are considered.
!! Only the turbulent fluxes of the conservative variables
-!! (Thetal and Rnp stored in PRx(:,:,:,1)) are computed.
-!! Note that the turbulent fluxes at the vertical
+!! (Thetal and Rnp stored in PRx(:,:,:,1)) are computed.
+!! Note that the turbulent fluxes at the vertical
!! boundaries are given either by the soil scheme for the surface one
-!! ( at the same instant as the others fluxes) and equal to 0 at the
-!! top of the model. The thermal production is computed by vertically
+!! ( at the same instant as the others fluxes) and equal to 0 at the
+!! top of the model. The thermal production is computed by vertically
!! averaging the turbulent flux and multiply this flux at the mass point by
!! a function ETHETA or EMOIST, which preform the transformation from the
-!! conservative variables to the virtual potential temperature.
-!!
+!! conservative variables to the virtual potential temperature.
+!!
!! In section 4, the variance of the statistical variable
-!! s indicating presence or not of condensation, is determined in function
+!! s indicating presence or not of condensation, is determined in function
!! of the turbulent moments of the conservative variables and its
-!! squarred root is stored in PSIGS. This information will be completed in
-!! the horizontal turbulence if the turbulence dimensionality is not
+!! squarred root is stored in PSIGS. This information will be completed in
+!! the horizontal turbulence if the turbulence dimensionality is not
!! equal to "1DIM".
!!
!! In section 5, the x component of the stress tensor is computed.
@@ -84,47 +84,47 @@ SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI, &
!! j" is also parallel to the surface and in the normal direction of
!! the maximum slope
!! k" is the normal to the surface
-!! In order to prevent numerical instability, the implicit scheme has
-!! been extended to the surface flux regarding to its dependence in
-!! function of U. The dependence in function of the other components
+!! In order to prevent numerical instability, the implicit scheme has
+!! been extended to the surface flux regarding to its dependence in
+!! function of U. The dependence in function of the other components
!! introduced by the different rotations is only explicit.
-!! The turbulent fluxes are used to compute the dynamic production of
+!! The turbulent fluxes are used to compute the dynamic production of
!! TKE. For the last TKE level ( located at PDZZ(:,:,IKB)/2 from the
-!! ground), an harmonic extrapolation from the dynamic production at
+!! ground), an harmonic extrapolation from the dynamic production at
!! PDZZ(:,:,IKB) is used to avoid an evaluation of the gradient of U
!! in the surface layer.
!!
!! In section 6, the same steps are repeated but for the y direction
-!! and in section 7, a diagnostic computation of the W variance is
+!! and in section 7, a diagnostic computation of the W variance is
!! performed.
!!
-!! In section 8, the turbulent fluxes for the scalar variables are
+!! In section 8, the turbulent fluxes for the scalar variables are
!! computed by the same way as the conservative thermodynamical variables
!!
-!!
+!!
!! EXTERNAL
!! --------
-!! GX_U_M, GY_V_M, GZ_W_M : cartesian gradient operators
+!! GX_U_M, GY_V_M, GZ_W_M : cartesian gradient operators
!! GX_U_UW,GY_V_VW (X,Y,Z) represent the direction of the gradient
-!! _(M,U,...)_ represent the localization of the
+!! _(M,U,...)_ represent the localization of the
!! field to be derivated
-!! _(M,UW,...) represent the localization of the
+!! _(M,UW,...) represent the localization of the
!! field derivated
-!!
+!!
!!
!! MXM,MXF,MYM,MYF,MZM,MZF
-!! : Shuman functions (mean operators)
+!! : Shuman functions (mean operators)
!! DXF,DYF,DZF,DZM
-!! : Shuman functions (difference operators)
-!!
-!! FUNCTIONs ETHETA and EMOIST :
+!! : Shuman functions (difference operators)
+!!
+!! FUNCTIONs ETHETA and EMOIST :
!! allows to compute:
!! - the coefficients for the turbulent correlation between
-!! any variable and the virtual potential temperature, of its
-!! correlations with the conservative potential temperature and
+!! any variable and the virtual potential temperature, of its
+!! correlations with the conservative potential temperature and
!! the humidity conservative variable:
!! ------- ------- -------
-!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
+!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
!!
!!
!! IMPLICIT ARGUMENTS
@@ -158,34 +158,34 @@ SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI, &
!! MODIFICATIONS
!! -------------
!! Original August 19, 1994
-!! Modifications: February 14, 1995 (J.Cuxart and J.Stein)
+!! Modifications: February 14, 1995 (J.Cuxart and J.Stein)
!! Doctorization and Optimization
-!! Modifications: March 21, 1995 (J.M. Carriere)
+!! Modifications: March 21, 1995 (J.M. Carriere)
!! Introduction of cloud water
-!! Modifications: June 14, 1995 (J.Cuxart and J. Stein)
+!! Modifications: June 14, 1995 (J.Cuxart and J. Stein)
!! Phi3 and Psi3 at w-point + bug in the all
-!! or nothing condens.
-!! Modifications: Sept 15, 1995 (J.Cuxart and J. Stein)
+!! or nothing condens.
+!! Modifications: Sept 15, 1995 (J.Cuxart and J. Stein)
!! Change the DP computation at the ground
-!! Modifications: October 10, 1995 (J.Cuxart and J. Stein)
+!! Modifications: October 10, 1995 (J.Cuxart and J. Stein)
!! Psi for scal var and LES tools
!! Modifications: November 10, 1995 (J. Stein)
-!! change the surface relations
+!! change the surface relations
!! Modifications: February 20, 1995 (J. Stein) optimization
-!! Modifications: May 21, 1996 (J. Stein)
-!! bug in the vertical flux of the V wind
+!! Modifications: May 21, 1996 (J. Stein)
+!! bug in the vertical flux of the V wind
!! component for explicit computation
-!! Modifications: May 21, 1996 (N. wood)
+!! Modifications: May 21, 1996 (N. wood)
!! modify the computation of the vertical
!! part or the surface tangential flux
!! Modifications: May 21, 1996 (P. Jabouille)
!! same modification in the Y direction
-!!
+!!
!! Modifications: Sept 17, 1996 (J. Stein) change the moist case by using
!! Pi instead of Piref + use Atheta and Amoist
!!
-!! Modifications: Nov 24, 1997 (V. Masson) removes the DO loops
-!! Modifications: Mar 31, 1998 (V. Masson) splits the routine TURB_VER_THERMO_FLUX
+!! Modifications: Nov 24, 1997 (V. Masson) removes the DO loops
+!! Modifications: Mar 31, 1998 (V. Masson) splits the routine TURB_VER_THERMO_FLUX
!! Modifications: Oct 18, 2000 (V. Masson) LES computations
!! Modifications: Dec 01, 2000 (V. Masson) conservation of energy from
!! surface flux in 1DIM case
@@ -193,12 +193,13 @@ SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI, &
!! Nov 06, 2002 (V. Masson) LES budgets
!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
!! change of YCOMMENT
-!! 2012-02 (Y. Seity) add possibility to run with reversed
+!! 2012-02 (Y. Seity) add possibility to run with reversed
!! vertical levels
!! Modifications July 2015 (Wim de Rooy) LHARAT switch
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Modifications June 2019 (Wim de Rooy) New set up cloud scheme
!!--------------------------------------------------------------------------
-!
+!
!* 0. DECLARATIONS
! ------------
!
@@ -230,7 +231,7 @@ IMPLICIT NONE
!
!
!
-INTEGER, INTENT(IN) :: KKA !near ground array index
+INTEGER, INTENT(IN) :: KKA !near ground array index
INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
INTEGER, INTENT(IN) :: KRR ! number of moist var.
@@ -252,31 +253,31 @@ REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSZW ! Director Cosinus of the
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry density * grid volum
REAL, DIMENSION(:,:,:), INTENT(IN) :: MFMOIST ! moist mass flux dual scheme
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! ref. state Virtual
- ! Potential Temperature
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! ref. state Virtual
+ ! Potential Temperature
!
REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHM,PSFRM ! surface fluxes at time
-! ! t - deltat
+! ! t - deltat
!
REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHP,PSFRP ! surface fluxes at time
-! ! t + deltat
+! ! t + deltat
!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM
! Vertical wind
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM
! potential temperature at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! Mixing ratios
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! Mixing ratios
! at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! Mixing ratios
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! Mixing ratios
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at time t
! In case LHARATU=TRUE, PLM already includes all stability corrections
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PLEPS ! dissipative length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: 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) :: PATHETA ! coefficients between
REAL, DIMENSION(:,:,:), INTENT(IN) :: PAMOIST ! s and Thetal and Rnp
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRCM ! normalized
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRCM ! normalized
! 2nd-order flux s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
REAL, DIMENSION(:,:,:), INTENT(IN) :: PBETA ! buoyancy coefficient
REAL, DIMENSION(:,:,:), INTENT(IN) :: PSQRT_TKE ! sqrt(e)
@@ -322,8 +323,8 @@ REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) :: &
PLMF, & ! estimate full level length scale from half levels (sub optimal)
PLEPSF ! estimate full level diss length scale from half levels (sub optimal)
-INTEGER :: IRESP ! Return code of FM routines
-INTEGER :: IGRID ! C-grid indicator in LFIFM file
+INTEGER :: IRESP ! Return code of FM routines
+INTEGER :: IGRID ! C-grid indicator in LFIFM file
INTEGER :: ILENCH ! Length of comment string in LFIFM file
INTEGER :: IKB,IKE ! I index values for the Beginning and End
INTEGER :: IKU ! array sizes
@@ -333,7 +334,7 @@ INTEGER :: I1,I2 ! For ZCOEFF allocation
CHARACTER (LEN=100) :: YCOMMENT ! comment string in LFIFM file
CHARACTER (LEN=16) :: YRECFM ! Name of the desired field in LFIFM file
REAL, DIMENSION(:,:,:),ALLOCATABLE :: ZCOEFF
- ! coefficients for the uncentred gradient
+ ! coefficients for the uncentred gradient
! computation near the ground
!
REAL :: ZTIME1, ZTIME2
@@ -363,7 +364,7 @@ ALLOCATE(ZCOEFF(SIZE(PDZZ,1),SIZE(PDZZ,2),I1:I2))
!
GUSERV = (KRR/=0)
!
-! compute the coefficients for the uncentred gradient computation near the
+! compute the coefficients for the uncentred gradient computation near the
! ground
ZCOEFF(:,:,IKB+2*KKL)= - PDZZ(:,:,IKB+KKL) / &
( (PDZZ(:,:,IKB+2*KKL)+PDZZ(:,:,IKB+KKL)) * PDZZ(:,:,IKB+2*KKL) )
@@ -374,15 +375,22 @@ ZCOEFF(:,:,IKB)= - (PDZZ(:,:,IKB+2*KKL)+2.*PDZZ(:,:,IKB+KKL)) / &
!
!
IF (LHARAT) THEN
-PLMF=MZF(PLM, KKA, KKU, KKL)
-PLEPSF=PLMF
-! function MZF produces -999 for level IKU (82 for 80 levels)
-! so put these to normal value as this level (82) is indeed calculated
-PLMF(:,:,IKU)=0.001
-PLEPSF(:,:,IKU)=0.001
-ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
+ PLMF=MZF(PLM, KKA, KKU, KKL)
+ !wc Part of the new statistical cloud scheme set up
+ IF (LSTATNW) THEN
+ PLEPSF=MZF(PLEPS, KKA, KKU, KKL)
+ ELSE
+ PLEPSF=PLMF
+ ENDIF
+ ! function MZF produces -999 for level IKU (82 for 80 levels)
+ ! so put these to normal value as this level (82) is indeed calculated
+ PLMF(:,:,IKU)=0.001
+ PLEPSF(:,:,IKU)=0.001
+ ! with energy cascade contribution 50MF term can be omitted
+ !ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) + 50*MFMOIST(:,:,:)
+ ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:))
ELSE
-ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
+ ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
ENDIF
!
@@ -409,15 +417,19 @@ END IF
! --------------------------------------------------------
!
!
-!* 4.2 <THl THl>
+!* 4.2 <THl THl>
!
! Compute the turbulent variance F and F' at time t-dt.
!
-IF (LHARAT) THEN
- ZF (:,:,:) = PLMF*PLEPSF*MZF(PDTH_DZ**2, KKA, KKU, KKL)
-ELSE
- ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(PPHI3*PDTH_DZ**2, KKA, KKU, KKL)
-ENDIF
+ IF (LHARAT) THEN
+ IF (LSTATNW) THEN
+ ZF (:,:,:) = XCTV*PLMF*PLEPSF*MZF(PDTH_DZ**2, KKA, KKU, KKL)
+ ELSE
+ ZF (:,:,:) = PLMF*PLEPSF*MZF(PDTH_DZ**2, KKA, KKU, KKL)
+ ENDIF
+ ELSE
+ ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(PPHI3*PDTH_DZ**2, KKA, KKU, KKL)
+ ENDIF
ZDFDDTDZ(:,:,:) = 0. ! this term, because of discretization, is treated separately
!
! Effect of 3rd order terms in temperature flux (at mass point)
@@ -473,41 +485,63 @@ ENDIF
!
! special case near the ground ( uncentred gradient )
IF (LHARAT) THEN
- ZFLXZ(:,:,IKB) = PLMF(:,:,IKB) &
- * PLEPSF(:,:,IKB) &
- *( PEXPL * &
- ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
- +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
- +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB ) )**2 &
- +PIMPL * &
- ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
- +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
- +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB ) )**2 &
- )
- ELSE
- ZFLXZ(:,:,IKB) = XCTV * PPHI3(:,:,IKB+KKL) * PLM(:,:,IKB) &
- * PLEPS(:,:,IKB) &
- *( PEXPL * &
- ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
- +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
- +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB ) )**2 &
- +PIMPL * &
- ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
- +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
- +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB ) )**2 &
- )
- ENDIF
+ IF (LSTATNW) THEN
+ ZFLXZ(:,:,IKB) = XCTV * PLMF(:,:,IKB) &
+ * PLEPSF(:,:,IKB) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB ) )**2 &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB ) )**2 &
+ )
+ ELSE
+ ZFLXZ(:,:,IKB) = PLMF(:,:,IKB) &
+ * PLEPSF(:,:,IKB) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB ) )**2 &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB ) )**2 &
+ )
+ ENDIF
+ ELSE
+ ZFLXZ(:,:,IKB) = XCTV * PPHI3(:,:,IKB+KKL) * PLM(:,:,IKB) &
+ * PLEPS(:,:,IKB) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB ) )**2 &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB ) )**2 &
+ )
+ ENDIF
!
- ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
!
- ZFLXZ = MAX(0., ZFLXZ)
+ IF (LSTATNW) THEN
+ !wc The variance from the budget eq should be multiplied by 2 here
+ ! thl'2=2*L*LEPS*(dthl/dz**2)
+ ZFLXZ = MAX(0., 2.*ZFLXZ)
+ ELSE
+ ZFLXZ = MAX(0., ZFLXZ)
+ ENDIF
!
IF (KRRL > 0) THEN
- PSIGS(:,:,:) = ZFLXZ(:,:,:) * PATHETA(:,:,:)**2
+ PSIGS(:,:,:) = ZFLXZ(:,:,:) * PATHETA(:,:,:)**2
+ ELSE
+ PSIGS(:,:,:) = 0.
END IF
!
!
- ! stores <THl THl>
+ ! stores <THl THl>
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
TZFIELD%CMNHNAME = 'THL_VVAR'
TZFIELD%CSTDNAME = ''
@@ -526,26 +560,30 @@ ENDIF
!
IF (LLES_CALL) THEN
CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID(ZFLXZ, X_LES_SUBGRID_Thl2 )
+ CALL LES_MEAN_SUBGRID(ZFLXZ, X_LES_SUBGRID_Thl2 )
CALL LES_MEAN_SUBGRID(MZF(PWM, KKA, KKU, KKL)*ZFLXZ, X_LES_RES_W_SBG_Thl2 )
- CALL LES_MEAN_SUBGRID(-2.*XCTD*PSQRT_TKE*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Thl2 )
- CALL LES_MEAN_SUBGRID(PETHETA*ZFLXZ, X_LES_SUBGRID_ThlThv )
- CALL LES_MEAN_SUBGRID(-XA3*PBETA*PETHETA*ZFLXZ, X_LES_SUBGRID_ThlPz, .TRUE. )
+ CALL LES_MEAN_SUBGRID(-2.*XCTD*PSQRT_TKE*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Thl2 )
+ CALL LES_MEAN_SUBGRID(PETHETA*ZFLXZ, X_LES_SUBGRID_ThlThv )
+ CALL LES_MEAN_SUBGRID(-XA3*PBETA*PETHETA*ZFLXZ, X_LES_SUBGRID_ThlPz, .TRUE. )
CALL SECOND_MNH(ZTIME2)
XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
END IF
!
IF ( KRR /= 0 ) THEN
!
-!* 4.3 <THl Rnp>
+!* 4.3 <THl Rnp>
!
!
! Compute the turbulent variance F and F' at time t-dt.
-IF (LHARAT) THEN
- ZF (:,:,:) = PLMF*PLEPSF*MZF(PDTH_DZ*PDR_DZ, KKA, KKU, KKL)
-ELSE
- ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(0.5*(PPHI3+PPSI3)*PDTH_DZ*PDR_DZ, KKA, KKU, KKL)
-ENDIF
+ IF (LHARAT) THEN
+ IF (LSTATNW) THEN
+ ZF (:,:,:) = XCTV*PLMF*PLEPSF*MZF(PDTH_DZ*PDR_DZ, KKA, KKU, KKL)
+ ELSE
+ ZF (:,:,:) = PLMF*PLEPSF*MZF(PDTH_DZ*PDR_DZ, KKA, KKU, KKL)
+ ENDIF
+ ELSE
+ ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(0.5*(PPHI3+PPSI3)*PDTH_DZ*PDR_DZ, KKA, KKU, KKL)
+ ENDIF
ZDFDDTDZ(:,:,:) = 0. ! this term, because of discretization, is treated separately
ZDFDDRDZ(:,:,:) = 0. ! this term, because of discretization, is treated separately
!
@@ -602,79 +640,119 @@ ENDIF
END IF
!
IF (LHARAT) THEN
- ZFLXZ(:,:,:) = ZF &
- + PIMPL * PLMF*PLEPSF*0.5 &
- * MZF(( 2. &
- ) *PDR_DZ *DZM(PTHLP - PTHLM, KKA, KKU, KKL ) / PDZZ &
- +( 2. &
- ) *PDTH_DZ *DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ &
- , KKA, KKU, KKL) &
- + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM(:,:,:), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
- + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
+ IF (LSTATNW) THEN
+ ZFLXZ(:,:,:) = ZF &
+ + PIMPL * XCTV*PLMF*PLEPSF*0.5 &
+ * MZF(( 2. &
+ ) *PDR_DZ *DZM(PTHLP - PTHLM, KKA, KKU, KKL ) / PDZZ &
+ +( 2. &
+ ) *PDTH_DZ *DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ &
+ , KKA, KKU, KKL) &
+ + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM(:,:,:), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
+ + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
+ ELSE
+ ZFLXZ(:,:,:) = ZF &
+ + PIMPL * PLMF*PLEPSF*0.5 &
+ * MZF(( 2. &
+ ) *PDR_DZ *DZM(PTHLP - PTHLM, KKA, KKU, KKL ) / PDZZ &
+ +( 2. &
+ ) *PDTH_DZ *DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ &
+ , KKA, KKU, KKL) &
+ + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM(:,:,:), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
+ + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
+ ENDIF
ELSE
- ZFLXZ(:,:,:) = ZF &
- + PIMPL * XCTV*PLM*PLEPS*0.5 &
- * MZF(( D_PHI3DTDZ_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV) & ! d(phi3*dthdz)/ddthdz term
- +D_PSI3DTDZ_O_DDTDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV) & ! d(psi3*dthdz)/ddthdz term
- ) *PDR_DZ *DZM(PTHLP - PTHLM, KKA, KKU, KKL) / PDZZ &
- +( D_PHI3DRDZ_O_DDRDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV) & ! d(phi3*drdz )/ddrdz term
- +D_PSI3DRDZ_O_DDRDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV) & ! d(psi3*drdz )/ddrdz term
- ) *PDTH_DZ *DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ &
- , KKA, KKU, KKL) &
- + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM(:,:,:), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
- + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
+ ZFLXZ(:,:,:) = ZF &
+ + PIMPL * XCTV*PLM*PLEPS*0.5 &
+ * MZF(( D_PHI3DTDZ_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV) & ! d(phi3*dthdz)/ddthdz term
+ +D_PSI3DTDZ_O_DDTDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV) & ! d(psi3*dthdz)/ddthdz term
+ ) *PDR_DZ *DZM(PTHLP - PTHLM, KKA, KKU, KKL) / PDZZ &
+ +( D_PHI3DRDZ_O_DDRDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV) & ! d(phi3*drdz )/ddrdz term
+ +D_PSI3DRDZ_O_DDRDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV) & ! d(psi3*drdz )/ddrdz term
+ ) *PDTH_DZ *DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ &
+ , KKA, KKU, KKL) &
+ + PIMPL * ZDFDDTDZ * MZF(DZM(PTHLP - PTHLM(:,:,:), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
+ + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
ENDIF
!
! special case near the ground ( uncentred gradient )
IF (LHARAT) THEN
- ZFLXZ(:,:,IKB) = &
- (1. ) &
- *( PEXPL * &
- ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
- +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
- +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB )) &
- *( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
- +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
- +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 )) &
- +PIMPL * &
- ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
- +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
- +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB )) &
- *( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL ) &
- +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
- +ZCOEFF(:,:,IKB )*PRP(:,:,IKB )) &
- )
- ELSE
- ZFLXZ(:,:,IKB) = &
- (XCHT1 * PPHI3(:,:,IKB+KKL) + XCHT2 * PPSI3(:,:,IKB+KKL)) &
- *( PEXPL * &
- ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
- +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
- +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB )) &
- *( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
- +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
- +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 )) &
- +PIMPL * &
- ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
- +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
- +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB )) &
- *( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL ) &
- +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
- +ZCOEFF(:,:,IKB )*PRP(:,:,IKB )) &
- )
+ IF (LSTATNW) THEN
+ ZFLXZ(:,:,IKB) = &
+ (XCHT1 + XCHT2 ) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB )) &
+ *( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
+ +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
+ +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 )) &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB )) &
+ *( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL ) &
+ +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PRP(:,:,IKB )) &
+ )
+ ELSE
+ ZFLXZ(:,:,IKB) = &
+ (1. ) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB )) &
+ *( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
+ +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
+ +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 )) &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB )) &
+ *( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL ) &
+ +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PRP(:,:,IKB )) &
+ )
+ ENDIF
+ ELSE
+ ZFLXZ(:,:,IKB) = &
+ (XCHT1 * PPHI3(:,:,IKB+KKL) + XCHT2 * PPSI3(:,:,IKB+KKL)) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLM(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLM(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB )) &
+ *( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
+ +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
+ +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 )) &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PTHLP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PTHLP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PTHLP(:,:,IKB )) &
+ *( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL ) &
+ +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PRP(:,:,IKB )) &
+ )
ENDIF
- !
- ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
!
- IF ( KRRL > 0 ) THEN
-!
-!
-! NB PATHETA is -b in Chaboureau Bechtold 2002 which explains the + sign here
-
- PSIGS(:,:,:) = PSIGS(:,:,:) + &
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ !
+ IF (LSTATNW) THEN
+ !wc The variance from the budget eq should be multiplied by 2 here
+ ! e.g. thl'2=2*L*LEPS*(cab)^-1 *(dthl/dz**2)
+ ZFLXZ = MIN(0., 2.*ZFLXZ)
+ ENDIF
+ IF ( KRRL > 0 ) THEN
+ IF (LSTATNW) THEN
+ !wc Part of the new statistical cloud scheme set up. Normal notation so - sign
+ PSIGS(:,:,:) = PSIGS(:,:,:) - &
2. * PATHETA(:,:,:) * PAMOIST(:,:,:) * ZFLXZ(:,:,:)
+ ELSE
+ ! NB PATHETA is -b in Chaboureau Bechtold 2002 which explains the + sign here
+ PSIGS(:,:,:) = PSIGS(:,:,:) + &
+ 2. * PATHETA(:,:,:) * PAMOIST(:,:,:) * ZFLXZ(:,:,:)
+ ENDIF
END IF
- ! stores <THl Rnp>
+ ! stores <THl Rnp>
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
TZFIELD%CMNHNAME = 'THLRCONS_VCOR'
TZFIELD%CSTDNAME = ''
@@ -693,27 +771,31 @@ ENDIF
!
IF (LLES_CALL) THEN
CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID(ZFLXZ, X_LES_SUBGRID_THlRt )
+ CALL LES_MEAN_SUBGRID(ZFLXZ, X_LES_SUBGRID_THlRt )
CALL LES_MEAN_SUBGRID(MZF(PWM, KKA, KKU, KKL)*ZFLXZ, X_LES_RES_W_SBG_ThlRt )
- CALL LES_MEAN_SUBGRID(-2.*XCTD*PSQRT_TKE*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_ThlRt )
- CALL LES_MEAN_SUBGRID(PETHETA*ZFLXZ, X_LES_SUBGRID_RtThv )
- CALL LES_MEAN_SUBGRID(-XA3*PBETA*PETHETA*ZFLXZ, X_LES_SUBGRID_RtPz, .TRUE. )
- CALL LES_MEAN_SUBGRID(PEMOIST*ZFLXZ, X_LES_SUBGRID_ThlThv , .TRUE. )
- CALL LES_MEAN_SUBGRID(-XA3*PBETA*PEMOIST*ZFLXZ, X_LES_SUBGRID_ThlPz, .TRUE. )
+ CALL LES_MEAN_SUBGRID(-2.*XCTD*PSQRT_TKE*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_ThlRt )
+ CALL LES_MEAN_SUBGRID(PETHETA*ZFLXZ, X_LES_SUBGRID_RtThv )
+ CALL LES_MEAN_SUBGRID(-XA3*PBETA*PETHETA*ZFLXZ, X_LES_SUBGRID_RtPz, .TRUE. )
+ CALL LES_MEAN_SUBGRID(PEMOIST*ZFLXZ, X_LES_SUBGRID_ThlThv , .TRUE. )
+ CALL LES_MEAN_SUBGRID(-XA3*PBETA*PEMOIST*ZFLXZ, X_LES_SUBGRID_ThlPz, .TRUE. )
CALL SECOND_MNH(ZTIME2)
XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
END IF
-!
+!
!
!* 4.4 <Rnp Rnp>
!
!
! Compute the turbulent variance F and F' at time t-dt.
-IF (LHARAT) THEN
- ZF (:,:,:) = PLMF*PLEPSF*MZF(PDR_DZ**2, KKA, KKU, KKL)
- ELSE
- ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(PPSI3*PDR_DZ**2, KKA, KKU, KKL)
-ENDIF
+ IF (LHARAT) THEN
+ IF (LSTATNW) THEN
+ ZF (:,:,:) = XCTV*PLMF*PLEPSF*MZF(PDR_DZ**2, KKA, KKU, KKL)
+ ELSE
+ ZF (:,:,:) = PLMF*PLEPSF*MZF(PDR_DZ**2, KKA, KKU, KKL)
+ ENDIF
+ ELSE
+ ZF (:,:,:) = XCTV*PLM*PLEPS*MZF(PPSI3*PDR_DZ**2, KKA, KKU, KKL)
+ ENDIF
ZDFDDRDZ(:,:,:) = 0. ! this term, because of discretization, is treated separately
!
! Effect of 3rd order terms in temperature flux (at mass point)
@@ -758,56 +840,83 @@ ENDIF
ZDFDDRDZ = ZDFDDRDZ + D_M3_R2_WTHR_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,&
& PD,PLEPS,PSQRT_TKE,PBLL_O_E,PETHETA,PDR_DZ) * PFTHR
END IF
-
+
END IF
!
- IF (LHARAT) THEN
- ZFLXZ(:,:,:) = ZF &
- + PIMPL * PLMF*PLEPSF &
- *MZF(2.*PDR_DZ* &
- DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
- + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
- ELSE
- ZFLXZ(:,:,:) = ZF &
- + PIMPL * XCTV*PLM*PLEPS &
- *MZF(D_PSI3DRDZ2_O_DDRDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,PDR_DZ,HTURBDIM,GUSERV) &
- *DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
- + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
- ENDIF
+ IF (LHARAT) THEN
+ IF (LSTATNW) THEN
+ ZFLXZ(:,:,:) = ZF &
+ + PIMPL * XCTV*PLMF*PLEPSF &
+ *MZF(2.*PDR_DZ* &
+ DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
+ + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
+ ELSE
+ ZFLXZ(:,:,:) = ZF &
+ + PIMPL * PLMF*PLEPSF &
+ *MZF(2.*PDR_DZ* &
+ DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
+ + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
+ ENDIF
+ ELSE
+ ZFLXZ(:,:,:) = ZF &
+ + PIMPL * XCTV*PLM*PLEPS &
+ *MZF(D_PSI3DRDZ2_O_DDRDZ(PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,PDR_DZ,HTURBDIM,GUSERV) &
+ *DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL) &
+ + PIMPL * ZDFDDRDZ * MZF(DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ, KKA, KKU, KKL)
+ ENDIF
!
! special case near the ground ( uncentred gradient )
- IF (LHARAT) THEN
- ZFLXZ(:,:,IKB) = PLMF(:,:,IKB) &
- * PLEPSF(:,:,IKB) &
- *( PEXPL * &
- ( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
- +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
- +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 ))**2 &
- +PIMPL * &
- ( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL) &
- +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
- +ZCOEFF(:,:,IKB )*PRP(:,:,IKB ))**2 &
- )
- ELSE
- ZFLXZ(:,:,IKB) = XCHV * PPSI3(:,:,IKB+KKL) * PLM(:,:,IKB) &
- * PLEPS(:,:,IKB) &
- *( PEXPL * &
- ( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
- +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
- +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 ))**2 &
- +PIMPL * &
- ( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL) &
- +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
- +ZCOEFF(:,:,IKB )*PRP(:,:,IKB ))**2 &
- )
- ENDIF
+ IF (LHARAT) THEN
+ IF (LSTATNW) THEN
+ ZFLXZ(:,:,IKB) = XCHV * PLMF(:,:,IKB) &
+ * PLEPSF(:,:,IKB) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
+ +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
+ +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 ))**2 &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PRP(:,:,IKB ))**2 &
+ )
+ ELSE
+ ZFLXZ(:,:,IKB) = PLMF(:,:,IKB) &
+ * PLEPSF(:,:,IKB) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
+ +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
+ +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 ))**2 &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PRP(:,:,IKB ))**2 &
+ )
+ ENDIF
+ ELSE
+ ZFLXZ(:,:,IKB) = XCHV * PPSI3(:,:,IKB+KKL) * PLM(:,:,IKB) &
+ * PLEPS(:,:,IKB) &
+ *( PEXPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PRM(:,:,IKB+2*KKL,1) &
+ +ZCOEFF(:,:,IKB+KKL )*PRM(:,:,IKB+KKL,1 ) &
+ +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1 ))**2 &
+ +PIMPL * &
+ ( ZCOEFF(:,:,IKB+2*KKL)*PRP(:,:,IKB+2*KKL) &
+ +ZCOEFF(:,:,IKB+KKL )*PRP(:,:,IKB+KKL ) &
+ +ZCOEFF(:,:,IKB )*PRP(:,:,IKB ))**2 &
+ )
+ ENDIF
!
- ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ IF (LSTATNW) THEN
+ !wc The variance from the budget eq should be multiplied by 2 here
+ ! thl'2=2*L*LEPS*(dthl/dz**2)
+ ZFLXZ = MAX(0., 2.*ZFLXZ)
+ ENDIF
!
IF ( KRRL > 0 ) THEN
PSIGS(:,:,:) = PSIGS(:,:,:) + PAMOIST(:,:,:) **2 * ZFLXZ(:,:,:)
END IF
- ! stores <Rnp Rnp>
+ ! stores <Rnp Rnp>
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
TZFIELD%CMNHNAME = 'RTOT_VVAR'
TZFIELD%CSTDNAME = ''
@@ -826,11 +935,11 @@ ENDIF
!
IF (LLES_CALL) THEN
CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID(ZFLXZ, X_LES_SUBGRID_Rt2 )
+ CALL LES_MEAN_SUBGRID(ZFLXZ, X_LES_SUBGRID_Rt2 )
CALL LES_MEAN_SUBGRID(MZF(PWM, KKA, KKU, KKL)*ZFLXZ, X_LES_RES_W_SBG_Rt2 )
- CALL LES_MEAN_SUBGRID(PEMOIST*ZFLXZ, X_LES_SUBGRID_RtThv , .TRUE. )
+ CALL LES_MEAN_SUBGRID(PEMOIST*ZFLXZ, X_LES_SUBGRID_RtThv , .TRUE. )
CALL LES_MEAN_SUBGRID(-XA3*PBETA*PEMOIST*ZFLXZ, X_LES_SUBGRID_RtPz, .TRUE. )
- CALL LES_MEAN_SUBGRID(-2.*XCTD*PSQRT_TKE*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Rt2 )
+ CALL LES_MEAN_SUBGRID(-2.*XCTD*PSQRT_TKE*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Rt2 )
CALL SECOND_MNH(ZTIME2)
XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
END IF
diff --git a/src/common/turb/mode_turb_ver_thermo_flux.F90 b/src/common/turb/mode_turb_ver_thermo_flux.F90
index e8966c17f42f22b2a492a7502316d08693e6ebda..4aac5c81ae2b9ec635b0f59654fc4c69924a2596 100644
--- a/src/common/turb/mode_turb_ver_thermo_flux.F90
+++ b/src/common/turb/mode_turb_ver_thermo_flux.F90
@@ -5,7 +5,7 @@
MODULE MODE_TURB_VER_THERMO_FLUX
IMPLICIT NONE
CONTAINS
-
+
SUBROUTINE TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI, &
OTURB_FLX,HTURBDIM,HTOM,OOCEAN, &
PIMPL,PEXPL, &
@@ -31,13 +31,13 @@ SUBROUTINE TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI, &
!! PURPOSE
!! -------
! The purpose of this routine is to compute the vertical turbulent
-! fluxes of the evolutive variables and give back the source
+! fluxes of the evolutive variables and give back the source
! terms to the main program. In the case of large horizontal meshes,
! the divergence of these vertical turbulent fluxes represent the whole
! effect of the turbulence but when the three-dimensionnal version of
! the turbulence scheme is activated (CTURBDIM="3DIM"), these divergences
-! are completed in the next routine TURB_HOR.
-! An arbitrary degree of implicitness has been implemented for the
+! are completed in the next routine TURB_HOR.
+! An arbitrary degree of implicitness has been implemented for the
! temporal treatment of these diffusion terms.
! The vertical boundary conditions are as follows:
! * at the bottom, the surface fluxes are prescribed at the same
@@ -45,8 +45,8 @@ SUBROUTINE TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI, &
! * at the top, the turbulent fluxes are set to 0.
! It should be noted that the condensation has been implicitely included
! in this turbulence scheme by using conservative variables and computing
-! the subgrid variance of a statistical variable s indicating the presence
-! or not of condensation in a given mesh.
+! the subgrid variance of a statistical variable s indicating the presence
+! or not of condensation in a given mesh.
!
!!** METHOD
!! ------
@@ -55,27 +55,27 @@ SUBROUTINE TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI, &
!! implicit scheme (a Crank-Nicholson type with coefficients different
!! than 0.5), which allows to vary the degree of implicitness of the
!! formulation.
-!! The different prognostic variables are treated one by one.
-!! The contributions of each turbulent fluxes are cumulated into the
-!! tendency PRvarS, and into the dynamic and thermal production of
+!! The different prognostic variables are treated one by one.
+!! The contributions of each turbulent fluxes are cumulated into the
+!! tendency PRvarS, and into the dynamic and thermal production of
!! TKE if necessary.
-!!
+!!
!! In section 2 and 3, the thermodynamical fields are considered.
!! Only the turbulent fluxes of the conservative variables
-!! (Thetal and Rnp stored in PRx(:,:,:,1)) are computed.
-!! Note that the turbulent fluxes at the vertical
+!! (Thetal and Rnp stored in PRx(:,:,:,1)) are computed.
+!! Note that the turbulent fluxes at the vertical
!! boundaries are given either by the soil scheme for the surface one
-!! ( at the same instant as the others fluxes) and equal to 0 at the
-!! top of the model. The thermal production is computed by vertically
+!! ( at the same instant as the others fluxes) and equal to 0 at the
+!! top of the model. The thermal production is computed by vertically
!! averaging the turbulent flux and multiply this flux at the mass point by
!! a function ETHETA or EMOIST, which preform the transformation from the
-!! conservative variables to the virtual potential temperature.
-!!
+!! conservative variables to the virtual potential temperature.
+!!
!! In section 4, the variance of the statistical variable
-!! s indicating presence or not of condensation, is determined in function
+!! s indicating presence or not of condensation, is determined in function
!! of the turbulent moments of the conservative variables and its
-!! squarred root is stored in PSIGS. This information will be completed in
-!! the horizontal turbulence if the turbulence dimensionality is not
+!! squarred root is stored in PSIGS. This information will be completed in
+!! the horizontal turbulence if the turbulence dimensionality is not
!! equal to "1DIM".
!!
!! In section 5, the x component of the stress tensor is computed.
@@ -86,53 +86,53 @@ SUBROUTINE TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI, &
!! j" is also parallel to the surface and in the normal direction of
!! the maximum slope
!! k" is the normal to the surface
-!! In order to prevent numerical instability, the implicit scheme has
-!! been extended to the surface flux regarding to its dependence in
-!! function of U. The dependence in function of the other components
+!! In order to prevent numerical instability, the implicit scheme has
+!! been extended to the surface flux regarding to its dependence in
+!! function of U. The dependence in function of the other components
!! introduced by the different rotations is only explicit.
-!! The turbulent fluxes are used to compute the dynamic production of
+!! The turbulent fluxes are used to compute the dynamic production of
!! TKE. For the last TKE level ( located at PDZZ(:,:,IKB)/2 from the
-!! ground), an harmonic extrapolation from the dynamic production at
+!! ground), an harmonic extrapolation from the dynamic production at
!! PDZZ(:,:,IKB) is used to avoid an evaluation of the gradient of U
!! in the surface layer.
!!
!! In section 6, the same steps are repeated but for the y direction
-!! and in section 7, a diagnostic computation of the W variance is
+!! and in section 7, a diagnostic computation of the W variance is
!! performed.
!!
-!! In section 8, the turbulent fluxes for the scalar variables are
+!! In section 8, the turbulent fluxes for the scalar variables are
!! computed by the same way as the conservative thermodynamical variables
!!
-!!
+!!
!! EXTERNAL
!! --------
-!! GX_U_M, GY_V_M, GZ_W_M : cartesian gradient operators
+!! GX_U_M, GY_V_M, GZ_W_M : cartesian gradient operators
!! GX_U_UW,GY_V_VW (X,Y,Z) represent the direction of the gradient
-!! _(M,U,...)_ represent the localization of the
+!! _(M,U,...)_ represent the localization of the
!! field to be derivated
-!! _(M,UW,...) represent the localization of the
+!! _(M,UW,...) represent the localization of the
!! field derivated
-!!
+!!
!!
!! MXM,MXF,MYM,MYF,MZM,MZF
-!! : Shuman functions (mean operators)
+!! : Shuman functions (mean operators)
!! DXF,DYF,DZF,DZM
-!! : Shuman functions (difference operators)
-!!
-!! SUBROUTINE TRIDIAG : to compute the split implicit evolution
+!! : Shuman functions (difference operators)
+!!
+!! SUBROUTINE TRIDIAG : to compute the splitted implicit evolution
!! of a variable located at a mass point
!!
-!! SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
+!! SUBROUTINE TRIDIAG_WIND: to compute the splitted implicit evolution
!! of a variable located at a wind point
!!
-!! FUNCTIONs ETHETA and EMOIST :
+!! FUNCTIONs ETHETA and EMOIST :
!! allows to compute:
!! - the coefficients for the turbulent correlation between
-!! any variable and the virtual potential temperature, of its
-!! correlations with the conservative potential temperature and
+!! any variable and the virtual potential temperature, of its
+!! correlations with the conservative potential temperature and
!! the humidity conservative variable:
!! ------- ------- -------
-!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
+!! A' Thv' = ETHETA A' Thl' + EMOIST A' Rnp'
!!
!!
!! IMPLICIT ARGUMENTS
@@ -166,34 +166,34 @@ SUBROUTINE TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI, &
!! MODIFICATIONS
!! -------------
!! Original August 19, 1994
-!! Modifications: February 14, 1995 (J.Cuxart and J.Stein)
+!! Modifications: February 14, 1995 (J.Cuxart and J.Stein)
!! Doctorization and Optimization
-!! Modifications: March 21, 1995 (J.M. Carriere)
+!! Modifications: March 21, 1995 (J.M. Carriere)
!! Introduction of cloud water
-!! Modifications: June 14, 1995 (J.Cuxart and J. Stein)
+!! Modifications: June 14, 1995 (J.Cuxart and J. Stein)
!! Phi3 and Psi3 at w-point + bug in the all
-!! or nothing condens.
-!! Modifications: Sept 15, 1995 (J.Cuxart and J. Stein)
+!! or nothing condens.
+!! Modifications: Sept 15, 1995 (J.Cuxart and J. Stein)
!! Change the DP computation at the ground
-!! Modifications: October 10, 1995 (J.Cuxart and J. Stein)
+!! Modifications: October 10, 1995 (J.Cuxart and J. Stein)
!! Psi for scal var and LES tools
!! Modifications: November 10, 1995 (J. Stein)
-!! change the surface relations
+!! change the surface relations
!! Modifications: February 20, 1995 (J. Stein) optimization
-!! Modifications: May 21, 1996 (J. Stein)
-!! bug in the vertical flux of the V wind
+!! Modifications: May 21, 1996 (J. Stein)
+!! bug in the vertical flux of the V wind
!! component for explicit computation
-!! Modifications: May 21, 1996 (N. wood)
+!! Modifications: May 21, 1996 (N. wood)
!! modify the computation of the vertical
!! part or the surface tangential flux
!! Modifications: May 21, 1996 (P. Jabouille)
!! same modification in the Y direction
-!!
+!!
!! Modifications: Sept 17, 1996 (J. Stein) change the moist case by using
!! Pi instead of Piref + use Atheta and Amoist
!!
-!! Modifications: Nov 24, 1997 (V. Masson) removes the DO loops
-!! Modifications: Mar 31, 1998 (V. Masson) splits the routine TURB_VER_THERMO_FLUX
+!! Modifications: Nov 24, 1997 (V. Masson) removes the DO loops
+!! Modifications: Mar 31, 1998 (V. Masson) splits the routine TURB_VER_THERMO_FLUX
!! Modifications: Oct 18, 2000 (V. Masson) LES computations
!! Modifications: Dec 01, 2000 (V. Masson) conservation of energy from
!! surface flux in 1DIM case
@@ -216,11 +216,13 @@ SUBROUTINE TURB_VER_THERMO_FLUX(KKA,KKU,KKL,KRR,KRRL,KRRI, &
!! applied to vertical fluxes of r_np and Thl
!! for implicit version of turbulence scheme
!! corrections and cleaning
+!! Modifications: June 2019 (Wim de Rooy) with energycascade, 50MF nog
+!! longer necessary
!! June 2020 (B. Vie) Patch preventing negative rc and ri in 2.3 and 3.3
!! JL Redelsperger : 03/2021: Ocean and Autocoupling O-A LES Cases
!! Sfc flux shape for LDEEPOC Case
!!--------------------------------------------------------------------------
-!
+!
!* 0. DECLARATIONS
! ------------
!
@@ -263,7 +265,7 @@ IMPLICIT NONE
!
!
!
-INTEGER, INTENT(IN) :: KKA !near ground array index
+INTEGER, INTENT(IN) :: KKA !near ground array index
INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
INTEGER, INTENT(IN) :: KRR ! number of moist var.
@@ -287,32 +289,32 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! altitudes
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry density * grid volum
REAL, DIMENSION(:,:,:), INTENT(IN) :: MFMOIST ! moist mass flux dual scheme
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! ref. state Virtual
- ! Potential Temperature
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! ref. state Virtual
+ ! Potential Temperature
!
REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHM,PSFRM ! surface fluxes at time
-! ! t - deltat
+! ! t - deltat
!
REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHP,PSFRP ! surface fluxes at time
-! ! t + deltat
+! ! t + deltat
!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM
! Vertical wind
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM
! potential temperature at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! Mixing ratios
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! Mixing ratios
! at t-Delta t
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! Mixing ratios
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVM ! Mixing ratios
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at time t
!
! In case LHARAT=TRUE, PLM already includes all stability corrections
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! Turb. mixing length
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PLEPS ! dissipative length
+REAL, DIMENSION(:,:,:), INTENT(IN) :: 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) :: PATHETA ! coefficients between
REAL, DIMENSION(:,:,:), INTENT(IN) :: PAMOIST ! s and Thetal and Rnp
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRCM ! normalized
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRCM ! normalized
! 2nd-order flux s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
REAL, DIMENSION(:,:,:), INTENT(IN) :: PFRAC_ICE ! ri fraction of rc+ri
REAL, DIMENSION(:,:,:), INTENT(IN) :: PBETA ! buoyancy coefficient
@@ -372,8 +374,8 @@ 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
-INTEGER :: JI, JJ ! loop indexes
+INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
+INTEGER :: JI, JJ ! loop indexes
!
INTEGER :: IIB,IJB ! Lower bounds of the physical
! sub-domain in x and y directions
@@ -404,23 +406,24 @@ TYPE(TFIELDDATA) :: TZFIELD
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB_VER_THERMO_FLUX',0,ZHOOK_HANDLE)
!
-! Size for a given proc & a given model
-IIU=SIZE(PTHLM,1)
+! Size for a given proc & a given model
+IIU=SIZE(PTHLM,1)
IJU=SIZE(PTHLM,2)
-IKT =SIZE(PTHLM,3)
-IKTE =IKT-JPVEXT_TURB
-IKTB =1+JPVEXT_TURB
+IKT =SIZE(PTHLM,3)
+IKTE =IKT-JPVEXT_TURB
+IKTB =1+JPVEXT_TURB
IKB=KKA+JPVEXT_TURB*KKL
IKE=KKU-JPVEXT_TURB*KKL
!
GUSERV = (KRR/=0)
!
-! compute the coefficients for the uncentred gradient computation near the
-! ground
+! compute the coefficients for the uncentred gradient computation near the ground
!
IF (LHARAT) THEN
-! LHARAT so TKE and length scales at half levels!
- ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) +50.*MFMOIST(:,:,:)
+ ! LHARAT so TKE and length scales at half levels!
+ !wc 50MF can be omitted with energy cascade included
+ !ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) +50.*MFMOIST(:,:,:)
+ ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:))
ELSE
ZKEFF(:,:,:) = MZM(PLM(:,:,:) * SQRT(PTKEM(:,:,:)), KKA, KKU, KKL)
ENDIF
@@ -454,8 +457,8 @@ IF (HTOM/='NONE') THEN
END IF
!----------------------------------------------------------------------------
!
-!* 2. SOURCES OF CONSERVATIVE POTENTIAL TEMPERATURE AND
-! PARTIAL THERMAL PRODUCTION
+!* 2. SOURCES OF CONSERVATIVE POTENTIAL TEMPERATURE AND
+! PARTIAL THERMAL PRODUCTION
! ---------------------------------------------------------------
!
!* 2.1 Splitted value for cons. potential temperature at t+deltat
@@ -527,13 +530,13 @@ IF (GFTHR) THEN
END IF
! compute interface flux
IF (LCOUPLES) THEN ! Autocoupling O-A LES
- IF (OOCEAN) THEN ! ocean model in coupled case
+ IF (OOCEAN) THEN ! ocean model in coupled case
ZF(:,:,IKE) = (XSSTFL_C(:,:,1)+XSSRFL_C(:,:,1)) &
*0.5* ( 1. + PRHODJ(:,:,KKU)/PRHODJ(:,:,IKE) )
ELSE ! atmosph model in coupled case
ZF(:,:,IKB) = XSSTFL_C(:,:,1) &
*0.5* ( 1. + PRHODJ(:,:,KKA)/PRHODJ(:,:,IKB) )
- ENDIF
+ ENDIF
!
ELSE ! No coupling O and A cases
! atmosp bottom
@@ -583,7 +586,7 @@ PRTHLS(:,:,:)= PRTHLS(:,:,:) + ZRWTHL(:,:,:)
!
!* 2.2 Partial Thermal Production
!
-! Conservative potential temperature flux :
+! Conservative potential temperature flux :
!
ZFLXZ(:,:,:) = ZF &
+ PIMPL * ZDFDDTDZ * DZM(PTHLP - PTHLM, KKA, KKU, KKL) / PDZZ
@@ -594,20 +597,20 @@ IF (LHGRAD) THEN
END WHERE
END IF
!
-ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
IF (OOCEAN) THEN
ZFLXZ(:,:,KKU) = ZFLXZ(:,:,IKE)
END IF
-!
+!
DO JK=IKTB+1,IKTE-1
PWTH(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+KKL))
END DO
!
-PWTH(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+KKL))
+PWTH(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+KKL))
!
IF (OOCEAN) THEN
PWTH(:,:,IKE)=0.5*(ZFLXZ(:,:,IKE)+ZFLXZ(:,:,IKE+KKL))
- PWTH(:,:,KKA)=0.
+ PWTH(:,:,KKA)=0.
PWTH(:,:,KKU)=ZFLXZ(:,:,KKU)
ELSE
PWTH(:,:,KKA)=0.5*(ZFLXZ(:,:,KKA)+ZFLXZ(:,:,KKA+KKL))
@@ -640,19 +643,19 @@ ELSE
ELSE
PTP(:,:,:)= PBETA * MZF( ZFLXZ,KKA, KKU, KKL )
END IF
-END IF
+END IF
!
! Buoyancy flux at flux points
-!
+!
PWTHV = MZM(PETHETA, KKA, KKU, KKL) * ZFLXZ
PWTHV(:,:,IKB) = PETHETA(:,:,IKB) * ZFLXZ(:,:,IKB)
!
IF (OOCEAN) THEN
- ! temperature contribution to Buy flux
+ ! temperature contribution to Buy flux
PWTHV(:,:,IKE) = PETHETA(:,:,IKE) * ZFLXZ(:,:,IKE)
END IF
!* 2.3 Partial vertical divergence of the < Rc w > flux
-! Correction for qc and qi negative in AROME
+! Correction for qc and qi negative in AROME
IF(CPROGRAM/='AROME ') THEN
IF ( KRRL >= 1 ) THEN
IF ( KRRI >= 1 ) THEN
@@ -670,16 +673,16 @@ IF(CPROGRAM/='AROME ') THEN
END IF
!
!* 2.4 Storage in LES configuration
-!
+!
IF (LLES_CALL) THEN
CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WThl )
+ CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WThl )
CALL LES_MEAN_SUBGRID(MZF(PWM*ZFLXZ, KKA, KKU, KKL), X_LES_RES_W_SBG_WThl )
CALL LES_MEAN_SUBGRID(GZ_W_M(PWM,PDZZ, KKA, KKU, KKL)*MZF(ZFLXZ, KKA, KKU, KKL),&
& X_LES_RES_ddxa_W_SBG_UaThl )
CALL LES_MEAN_SUBGRID(MZF(PDTH_DZ*ZFLXZ, KKA, KKU, KKL), X_LES_RES_ddxa_Thl_SBG_UaThl )
- CALL LES_MEAN_SUBGRID(-XCTP*PSQRT_TKE/PLM*MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_ThlPz )
- CALL LES_MEAN_SUBGRID(MZF(MZM(PETHETA, KKA, KKU, KKL)*ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WThv )
+ CALL LES_MEAN_SUBGRID(-XCTP*PSQRT_TKE/PLM*MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_ThlPz )
+ CALL LES_MEAN_SUBGRID(MZF(MZM(PETHETA, KKA, KKU, KKL)*ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WThv )
IF (KRR>=1) THEN
CALL LES_MEAN_SUBGRID(MZF(PDR_DZ*ZFLXZ, KKA, KKU, KKL), X_LES_RES_ddxa_Rt_SBG_UaThl )
END IF
@@ -690,21 +693,21 @@ IF (LLES_CALL) THEN
ZA(:,:,IKB) = XCSHF*PPHI3(:,:,IKB)*ZKEFF(:,:,IKB)
ZA = MZF(ZA, KKA, KKU, KKL)
ZA = MIN(MAX(ZA,-1000.),1000.)
- CALL LES_MEAN_SUBGRID( ZA, X_LES_SUBGRID_Kh )
+ CALL LES_MEAN_SUBGRID( ZA, X_LES_SUBGRID_Kh )
!
CALL SECOND_MNH(ZTIME2)
XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
END IF
!
!* 2.5 New boundary layer depth for TOMs
-!
+!
IF (HTOM=='TM06') CALL TM06_H(IKB,IKTB,IKTE,PTSTEP,PZZ,ZFLXZ,PBL_DEPTH)
!
!----------------------------------------------------------------------------
!
!
-!* 3. SOURCES OF CONSERVATIVE AND CLOUD MIXING RATIO AND
-! COMPLETE THERMAL PRODUCTION
+!* 3. SOURCES OF CONSERVATIVE AND CLOUD MIXING RATIO AND
+! COMPLETE THERMAL PRODUCTION
! ------------------------------------------------------
!
!* 3.1 Splitted value for cons. mixing ratio at t+deltat
@@ -755,7 +758,7 @@ IF (KRR /= 0) THEN
IF (GFWTH) THEN
ZF = ZF + M3_WR_W2TH(KKA,KKU,KKL,PD,ZKEFF,&
& PTKEM,PBLL_O_E,PETHETA,PDR_DZ) * PFWTH
- ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_W2TH_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,&
+ ZDFDDRDZ = ZDFDDRDZ + D_M3_WR_W2TH_O_DDRDZ(KKA,KKU,KKL,PREDR1,PREDTH1,&
& PD,ZKEFF,PTKEM,PBLL_O_E,PETHETA) * PFWTH
END IF
!
@@ -841,7 +844,7 @@ IF (KRR /= 0) THEN
! cons. mixing ratio flux :
!
ZFLXZ(:,:,:) = ZF &
- + PIMPL * ZDFDDRDZ * DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ
+ + PIMPL * ZDFDDRDZ * DZM(PRP - PRM(:,:,:,1), KKA, KKU, KKL) / PDZZ
!
! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
IF (LHGRAD) THEN
@@ -850,7 +853,7 @@ IF (KRR /= 0) THEN
END WHERE
END IF
!
- ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
!
DO JK=IKTB+1,IKTE-1
PWRC(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+KKL))
@@ -886,12 +889,12 @@ IF (KRR /= 0) THEN
END IF
!
! Buoyancy flux at flux points
- !
+ !
PWTHV = PWTHV + MZM(PEMOIST, KKA, KKU, KKL) * ZFLXZ
PWTHV(:,:,IKB) = PWTHV(:,:,IKB) + PEMOIST(:,:,IKB) * ZFLXZ(:,:,IKB)
IF (OOCEAN) THEN
PWTHV(:,:,IKE) = PWTHV(:,:,IKE) + PEMOIST(:,:,IKE)* ZFLXZ(:,:,IKE)
- END IF
+ END IF
!
!* 3.3 Complete vertical divergence of the < Rc w > flux
! Correction of qc and qi negative for AROME
@@ -912,17 +915,17 @@ IF(CPROGRAM/='AROME ') THEN
END IF
!
!* 3.4 Storage in LES configuration
-!
+!
IF (LLES_CALL) THEN
CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WRt )
+ CALL LES_MEAN_SUBGRID(MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WRt )
CALL LES_MEAN_SUBGRID(MZF(PWM*ZFLXZ, KKA, KKU, KKL), X_LES_RES_W_SBG_WRt )
CALL LES_MEAN_SUBGRID(GZ_W_M(PWM,PDZZ, KKA, KKU, KKL)*MZF(ZFLXZ, KKA, KKU, KKL),&
& X_LES_RES_ddxa_W_SBG_UaRt )
CALL LES_MEAN_SUBGRID(MZF(PDTH_DZ*ZFLXZ, KKA, KKU, KKL), X_LES_RES_ddxa_Thl_SBG_UaRt )
CALL LES_MEAN_SUBGRID(MZF(PDR_DZ*ZFLXZ, KKA, KKU, KKL), X_LES_RES_ddxa_Rt_SBG_UaRt )
- CALL LES_MEAN_SUBGRID(MZF(MZM(PEMOIST, KKA, KKU, KKL)*ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WThv , .TRUE. )
- CALL LES_MEAN_SUBGRID(-XCTP*PSQRT_TKE/PLM*MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_RtPz )
+ CALL LES_MEAN_SUBGRID(MZF(MZM(PEMOIST, KKA, KKU, KKL)*ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WThv , .TRUE. )
+ CALL LES_MEAN_SUBGRID(-XCTP*PSQRT_TKE/PLM*MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_RtPz )
CALL SECOND_MNH(ZTIME2)
XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
END IF
@@ -936,28 +939,28 @@ END IF
! -------------------------------
!
!
-!* 4.1 <w Rc>
+!* 4.1 <w Rc>
!
IF ( ((OTURB_FLX .AND. TPFILE%LOPENED) .OR. LLES_CALL) .AND. (KRRL > 0) ) THEN
-!
-! recover the Conservative potential temperature flux :
+!
+! recover the Conservative potential temperature flux :
! With LHARAT is true tke and length scales at half levels
! yet modify to use length scale and tke at half levels from vdfexcuhl
IF (LHARAT) THEN
ZA(:,:,:) = DZM(PIMPL * PTHLP + PEXPL * PTHLM, KKA, KKU, KKL) / PDZZ * &
- (-PLM*PSQRT_TKE)
+ (-PLM*PSQRT_TKE)
ELSE
ZA(:,:,:) = DZM(PIMPL * PTHLP + PEXPL * PTHLM, KKA, KKU, KKL) / PDZZ * &
- (-PPHI3*MZM(PLM*PSQRT_TKE, KKA, KKU, KKL)) * XCSHF
+ (-PPHI3*MZM(PLM*PSQRT_TKE, KKA, KKU, KKL)) * XCSHF
ENDIF
ZA(:,:,IKB) = ( PIMPL*PSFTHP(:,:) + PEXPL*PSFTHM(:,:) ) &
* PDIRCOSZW(:,:)
- !
+ !
! compute <w Rc>
ZFLXZ(:,:,:) = MZM(PAMOIST * 2.* PSRCM, KKA, KKU, KKL) * ZFLXZ(:,:,:) + &
MZM(PATHETA * 2.* PSRCM, KKA, KKU, KKL) * ZA(:,:,:)
- ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
- !
+ ZFLXZ(:,:,KKA) = ZFLXZ(:,:,IKB)
+ !
! store the liquid water mixing ratio vertical flux
IF ( OTURB_FLX .AND. TPFILE%LOPENED ) THEN
TZFIELD%CMNHNAME = 'RCW_FLX'
@@ -972,12 +975,12 @@ IF ( ((OTURB_FLX .AND. TPFILE%LOPENED) .OR. LLES_CALL) .AND. (KRRL > 0) ) THEN
TZFIELD%LTIMEDEP = .TRUE.
CALL IO_Field_write(TPFILE,TZFIELD,ZFLXZ)
END IF
- !
+ !
! and we store in LES configuration this subgrid flux <w'rc'>
!
IF (LLES_CALL) THEN
CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID( MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WRc )
+ CALL LES_MEAN_SUBGRID( MZF(ZFLXZ, KKA, KKU, KKL), X_LES_SUBGRID_WRc )
CALL SECOND_MNH(ZTIME2)
XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
END IF
diff --git a/src/common/turb/turb.F90 b/src/common/turb/turb.F90
index f66202783fd688fb91f30dedf586653e4e546d11..63fd923053c57e7c79e2bcdd40e1ff242c83d09e 100644
--- a/src/common/turb/turb.F90
+++ b/src/common/turb/turb.F90
@@ -24,26 +24,26 @@
& PEDR,PLEM,PRTKEMS,PTPMF, &
& PDRUS_TURB,PDRVS_TURB, &
& PDRTHLS_TURB,PDRRTS_TURB,PDRSVS_TURB,PTR,PDISS, &
- & PCURRENT_TKE_DISS )
+ & PCURRENT_TKE_DISS )
! #################################################################
!
!
!!**** *TURB* - computes the turbulent source terms for the prognostic
-!! variables.
+!! variables.
!!
!! PURPOSE
!! -------
-!!**** The purpose of this routine is to compute the source terms in
-!! the evolution equations due to the turbulent mixing.
+!!**** The purpose of this routine is to compute the source terms in
+!! the evolution equations due to the turbulent mixing.
!! The source term is computed as the divergence of the turbulent fluxes.
!! The cartesian fluxes are obtained by a one and a half order closure, based
-!! on a prognostic equation for the Turbulence Kinetic Energy( TKE ). The
-!! system is closed by prescribing a turbulent mixing length. Different
-!! choices are available for this length.
+!! on a prognostic equation for the Turbulence Kinetic Energy( TKE ). The
+!! system is closed by prescribing a turbulent mixing length. Different
+!! choices are available for this length.
!
!!** METHOD
!! ------
-!!
+!!
!! The dimensionality of the turbulence parameterization can be chosen by
!! means of the parameter HTURBDIM:
!! * HTURBDIM='1DIM' the parameterization is 1D but can be used in
@@ -51,7 +51,7 @@
!! turbulent fluxes are taken into account.
!! * HTURBDIM='3DIM' the parameterization is fully 2D or 3D depending
!! on the model dimensionality. Of course, it does not make any sense to
-!! activate this option with a 1D model.
+!! activate this option with a 1D model.
!!
!! The following steps are made:
!! 1- Preliminary computations.
@@ -62,39 +62,39 @@
!! original level of an air particule having an initial internal
!! energy equal to its TKE and stopped by the buoyancy forces.
!! The discrete formulation is second order accurate.
-!! * HTURBLEN='DELT' the mixing length is given by the mesh size
-!! depending on the model dimensionality, this length is limited
+!! * HTURBLEN='DELT' the mixing length is given by the mesh size
+!! depending on the model dimensionality, this length is limited
!! with the ground distance.
-!! * HTURBLEN='DEAR' the mixing length is given by the mesh size
-!! depending on the model dimensionality, this length is limited
+!! * HTURBLEN='DEAR' the mixing length is given by the mesh size
+!! depending on the model dimensionality, this length is limited
!! with the ground distance and also by the Deardorff mixing length
!! pertinent in the stable cases.
-!! * HTURBLEN='KEPS' the mixing length is deduced from the TKE
+!! * HTURBLEN='KEPS' the mixing length is deduced from the TKE
!! dissipation, which becomes a prognostic variable of the model (
-!! Duynkerke formulation).
+!! Duynkerke formulation).
!! 3'- The cloud mixing length is computed according to HTURBLEN_CLOUD
!! and emphasized following the CEI index
!! 4- The conservative variables are computed along with Lv/Cp.
!! 5- The turbulent Prandtl numbers are computed from the resolved fields
-!! and TKE
+!! and TKE
!! 6- The sources associated to the vertical turbulent fluxes are computed
-!! with a temporal scheme allowing a degree of implicitness given by
+!! with a temporal scheme allowing a degree of implicitness given by
!! PIMPL, varying from PIMPL=0. ( purely explicit scheme) to PIMPL=1.
!! ( purely implicit scheme)
!! The sources associated to the horizontal fluxes are computed with a
!! purely explicit temporal scheme. These sources are only computed when
!! the turbulence parameterization is 2D or 3D( HTURBDIM='3DIM' ).
-!! 7- The sources for TKE are computed, along with the dissipation of TKE
+!! 7- The sources for TKE are computed, along with the dissipation of TKE
!! if HTURBLEN='KEPS'.
-!! 8- Some turbulence-related quantities are stored in the synchronous
+!! 8- Some turbulence-related quantities are stored in the synchronous
!! FM-file.
-!! 9- The non-conservative variables are retrieved.
-!!
-!!
+!! 9- The non-conservative variables are retrieved.
+!!
+!!
!! The saving of the fields in the synchronous FM-file is controlled by:
!! * OTURB_FLX => saves all the turbulent fluxes and correlations
!! * OTURB_DIAG=> saves the turbulent Prandtl and Schmidt numbers, the
-!! source terms of TKE and dissipation of TKE
+!! source terms of TKE and dissipation of TKE
!!
!! EXTERNAL
!! --------
@@ -121,24 +121,24 @@
!!
!! MODD_CTURB : contains turbulence scheme constants
!! XCMFS,XCED to compute the dissipation mixing length
-!! XTKEMIN minimum values for the TKE
-!! XLINI,XLINF to compute Bougeault-Lacarrere mixing
+!! XTKEMIN minimum values for the TKE
+!! XLINI,XLINF to compute Bougeault-Lacarrere mixing
!! length
!! Module MODD_BUDGET:
-!! NBUMOD
-!! CBUTYPE
-!! LBU_RU
-!! LBU_RV
-!! LBU_RW
-!! LBU_RTH
-!! LBU_RSV1
-!! LBU_RRV
-!! LBU_RRC
-!! LBU_RRR
-!! LBU_RRI
-!! LBU_RRS
-!! LBU_RRG
-!! LBU_RRH
+!! NBUMOD
+!! CBUTYPE
+!! LBU_RU
+!! LBU_RV
+!! LBU_RW
+!! LBU_RTH
+!! LBU_RSV1
+!! LBU_RRV
+!! LBU_RRC
+!! LBU_RRR
+!! LBU_RRI
+!! LBU_RRS
+!! LBU_RRG
+!! LBU_RRH
!!
!! REFERENCE
!! ---------
@@ -152,11 +152,11 @@
!! MODIFICATIONS
!! -------------
!! Original 05/10/94
-!! Modifications: Feb 14, 1995 (J.Cuxart and J.Stein)
+!! Modifications: Feb 14, 1995 (J.Cuxart and J.Stein)
!! Doctorization and Optimization
-!! Modifications: March 21, 1995 (J.M. Carriere)
+!! Modifications: March 21, 1995 (J.M. Carriere)
!! Introduction of cloud water
-!! Modifications: June 1, 1995 (J.Cuxart )
+!! Modifications: June 1, 1995 (J.Cuxart )
!! take min(Kz,delta)
!! Modifications: June 1, 1995 (J.Stein J.Cuxart)
!! remove unnecessary arrays and change Prandtl
@@ -164,18 +164,18 @@
!! Modifications: July 20, 1995 (J.Stein) remove MODI_ground_ocean +
!! TZDTCUR + MODD_TIME because they are not used
!! change RW in RNP for the outputs
-!! Modifications: August 21, 1995 (Ph. Bougeault)
+!! Modifications: August 21, 1995 (Ph. Bougeault)
!! take min(K(z-zsol),delta)
!! Modifications: Sept 14, 1995 (Ph Bougeault, J. Cuxart)
-!! second order BL89 mixing length computations + add Deardorff length
+!! second order BL89 mixing length computations + add Deardorff length
!! in the Delta case for stable cases
!! Modifications: Sept 19, 1995 (J. Stein, J. Cuxart)
!! define a DEAR case for the mixing length, add MODI_BUDGET and change
!! some BUDGET calls, add LES tools
!! Modifications: Oct 16, 1995 (J. Stein) change the budget calls
-!! Modifications: Feb 28, 1996 (J. Stein) optimization +
+!! Modifications: Feb 28, 1996 (J. Stein) optimization +
!! remove min(K(z-zsol),delta)+
-!! bug in the tangential fluxes
+!! bug in the tangential fluxes
!! Modifications: Oct 16, 1996 (J. Stein) change the subgrid condensation
!! scheme + temporal discretization
!! Modifications: Dec 19, 1996 (J.-P. Pinty) update the budget calls
@@ -198,7 +198,7 @@
!! Feb 20, 2003 (J.-P. Pinty) Add reversible ice processes
!! May,26 2004 (P Jabouille) coef for computing dissipative heating
!! Sept 2004 (M.Tomasini) Cloud Mixing length modification
-!! following the instability
+!! following the instability
!! criterium CEI calculated in modeln
!! May 2006 Remove KEPS
!! Sept.2006 (I.Sandu): Modification of the stability criterion for
@@ -220,6 +220,7 @@
! P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
! Q. Rodier 01/2018: introduction of RM17
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!! June 2019 (Wim de Rooy) update statistical cloud scheme
! P. Wautelet 02/2020: use the new data structures and subroutines for budgets
! B. Vie 03/2020: LIMA negativity checks after turbulence, advection and microphysics budgets
! P. Wautelet 11/06/2020: bugfix: correct PRSVS array indices
@@ -283,10 +284,10 @@ IMPLICIT NONE
!
!
!
-INTEGER, INTENT(IN) :: KKA !near ground array index
+INTEGER, INTENT(IN) :: KKA !near ground array index
INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
-INTEGER, INTENT(IN) :: KMI ! model index number
+INTEGER, INTENT(IN) :: KMI ! model index number
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.
@@ -297,7 +298,7 @@ LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
! turbulent fluxes in the syncronous FM-file
LOGICAL, INTENT(IN) :: OTURB_DIAG ! switch to write some
! diagnostic fields in the syncronous FM-file
-LOGICAL, INTENT(IN) :: OSUBG_COND ! switch for SUBGrid
+LOGICAL, INTENT(IN) :: OSUBG_COND ! switch for SUBGrid
! CONDensation
LOGICAL, INTENT(IN) :: ORMC01 ! switch for RMC01 lengths in SBL
LOGICAL, INTENT(IN) :: OOCEAN ! switch for Ocean model version
@@ -309,12 +310,12 @@ CHARACTER(LEN=4), INTENT(IN) :: HTOM ! kind of Third Order Mome
CHARACTER(LEN=4), INTENT(IN) :: HTURBLEN_CL ! kind of cloud mixing length
REAL, INTENT(IN) :: PIMPL ! degree of implicitness
CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! Kind of microphysical scheme
-REAL, INTENT(IN) :: PTSTEP ! timestep
+REAL, INTENT(IN) :: PTSTEP ! timestep
TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ,PDZX,PDZY
! metric coefficients
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! physical distance
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! physical distance
! between 2 succesive grid points along the K direction
REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSXW, PDIRCOSYW, PDIRCOSZW
! Director Cosinus along x, y and z directions at surface w-point
@@ -328,11 +329,11 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! Virtual Potential
! Temperature of the reference state
!
REAL, DIMENSION(:,:), INTENT(IN) :: PSFTH,PSFRV, &
-! normal surface fluxes of theta and Rv
+! normal surface fluxes of theta and Rv
PSFU,PSFV
! normal surface fluxes of (u,v) parallel to the orography
REAL, DIMENSION(:,:,:), INTENT(IN) :: PSFSV
-! normal surface fluxes of Scalar var.
+! normal surface fluxes of Scalar var.
!
! prognostic variables at t- deltat
REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! Pressure at time t
@@ -346,7 +347,7 @@ REAL, DIMENSION(:,:), INTENT(INOUT) :: PSBL_DEPTH ! SBL depth for RMC01
!
! variables for cloud mixing length
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCEI ! Cloud Entrainment instability
- ! index to emphasize localy
+ ! index to emphasize localy
! turbulent fluxes
REAL, INTENT(IN) :: PCEI_MIN ! minimum threshold for the instability index CEI
REAL, INTENT(IN) :: PCEI_MAX ! maximum threshold for the instability index CEI
@@ -354,10 +355,10 @@ REAL, INTENT(IN) :: PCOEF_AMPL_SAT ! saturation of the amplification coeff
!
! thermodynamical variables which are transformed in conservative var.
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHLT ! conservative pot. temp.
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRT ! water var. where
- ! PRT(:,:,:,1) is the conservative mixing ratio
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRT ! water var. where
+ ! PRT(:,:,:,1) is the conservative mixing ratio
!
-! sources of momentum, conservative potential temperature, Turb. Kin. Energy,
+! sources of momentum, conservative potential temperature, Turb. Kin. Energy,
! TKE dissipation
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS,PRVS,PRWS,PRTHLS,PRTKES
! Source terms for all water kinds, PRRS(:,:,:,1) is used for the conservative
@@ -366,7 +367,7 @@ REAL, DIMENSION(:,:,:), INTENT(IN),OPTIONAL :: PRTKEMS
REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS
! Source terms for all passive scalar variables
REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRSVS
-! Sigma_s at time t+1 : square root of the variance of the deviation to the
+! Sigma_s at time t+1 : square root of the variance of the deviation to the
! saturation
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSIGS
REAL, DIMENSION(:,:,:), INTENT(OUT),OPTIONAL :: PDRUS_TURB ! evolution of rhoJ*U by turbulence only
@@ -374,7 +375,7 @@ REAL, DIMENSION(:,:,:), INTENT(OUT),OPTIONAL :: PDRVS_TURB ! evolution of
REAL, DIMENSION(:,:,:), INTENT(OUT),OPTIONAL :: PDRTHLS_TURB ! evolution of rhoJ*thl by turbulence only
REAL, DIMENSION(:,:,:), INTENT(OUT),OPTIONAL :: PDRRTS_TURB ! evolution of rhoJ*rt by turbulence only
REAL, DIMENSION(:,:,:,:), INTENT(OUT),OPTIONAL :: PDRSVS_TURB ! evolution of rhoJ*Sv by turbulence only
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PFLXZTHVMF
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PFLXZTHVMF
! MF contribution for vert. turb. transport
! used in the buoy. prod. of TKE
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PWTH ! heat flux
@@ -412,7 +413,7 @@ REAL, DIMENSION(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)) :: &
ZLOCPEXNM, & ! Lv/Cp/EXNREF at t-1
ZLM,ZLMW, & ! Turbulent mixing length (+ work array)
ZLEPS, & ! Dissipative length
- ZTRH, & !
+ ZTRH, & !
ZATHETA,ZAMOIST, & ! coefficients for s = f (Thetal,Rnp)
ZCOEF_DISS, & ! 1/(Cph*Exner) for dissipative heating
ZFRAC_ICE, & ! ri fraction of rc+ri
@@ -420,15 +421,15 @@ REAL, DIMENSION(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)) :: &
ZFWTH,ZFWR,ZFTH2,ZFR2,ZFTHR,& ! opposite of verticale derivate of 3rd order moments
ZTHLM,ZRTKEMS ! initial potential temp; TKE advective source
REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),SIZE(PRT,4)) :: &
- ZRM ! initial mixing ratio
+ ZRM ! initial mixing ratio
REAL, DIMENSION(SIZE(PTHLT,1),SIZE(PTHLT,2)) :: ZTAU11M,ZTAU12M, &
ZTAU22M,ZTAU33M, &
! tangential surface fluxes in the axes following the orography
ZUSLOPE,ZVSLOPE, &
- ! wind components at the first mass level parallel
- ! to the orography
+ ! wind components at the first mass level parallel
+ ! to the orography
ZCDUEFF, &
- ! - Cd*||u|| where ||u|| is the module of the wind tangential to
+ ! - Cd*||u|| where ||u|| is the module of the wind tangential to
! orography (ZUSLOPE,ZVSLOPE) at the surface.
ZUSTAR, ZLMO, &
ZRVM, ZSFRV
@@ -436,7 +437,7 @@ REAL, DIMENSION(SIZE(PTHLT,1),SIZE(PTHLT,2)) :: ZTAU11M,ZTAU12M, &
!
! Virtual Potential Temp. used
! in the Deardorff mixing length computation
-REAL, DIMENSION(:,:,:), ALLOCATABLE :: &
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: &
ZLVOCPEXNM,ZLSOCPEXNM, & ! Lv/Cp/EXNREF and Ls/Cp/EXNREF at t-1
ZATHETA_ICE,ZAMOIST_ICE ! coefficients for s = f (Thetal,Rnp)
!
@@ -446,12 +447,12 @@ REAL :: ZRVORD ! RV/RD
INTEGER :: IKB,IKE ! index value for the
! Beginning and the End of the physical domain for the mass points
INTEGER :: IKT ! array size in k direction
-INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
+INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
INTEGER :: JRR,JK,JSV ! loop counters
INTEGER :: JI,JJ ! loop counters
REAL :: ZL0 ! Max. Mixing Length in Blakadar formula
-REAL :: ZALPHA ! work coefficient :
- ! - proportionnality constant between Dz/2 and
+REAL :: ZALPHA ! work coefficient :
+ ! - proportionnality constant between Dz/2 and
! ! BL89 mixing length near the surface
!
REAL :: ZTIME1, ZTIME2
@@ -461,7 +462,7 @@ TYPE(TFIELDDATA) :: TZFIELD
!* 1.PRELIMINARIES
! -------------
!
-!* 1.1 Set the internal domains, ZEXPL
+!* 1.1 Set the internal domains, ZEXPL
!
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
@@ -474,8 +475,8 @@ IF (LHARAT .AND. LLES_CALL) THEN
CALL ABOR1('LHARATU not implemented for option LLES_CALL')
ENDIF
-IKT=SIZE(PTHLT,3)
-IKTB=1+JPVEXT_TURB
+IKT=SIZE(PTHLT,3)
+IKTB=1+JPVEXT_TURB
IKTE=IKT-JPVEXT_TURB
IKB=KKA+JPVEXT_TURB*KKL
IKE=KKU-JPVEXT_TURB*KKL
@@ -502,11 +503,11 @@ END IF
ZCP(:,:,:)=XCPD
!
IF (KRR > 0) ZCP(:,:,:) = ZCP(:,:,:) + XCPV * PRT(:,:,:,1)
-DO JRR = 2,1+KRRL ! loop on the liquid components
+DO JRR = 2,1+KRRL ! loop on the liquid components
ZCP(:,:,:) = ZCP(:,:,:) + XCL * PRT(:,:,:,JRR)
END DO
!
-DO JRR = 2+KRRL,1+KRRL+KRRI ! loop on the solid components
+DO JRR = 2+KRRL,1+KRRL+KRRI ! loop on the solid components
ZCP(:,:,:) = ZCP(:,:,:) + XCI * PRT(:,:,:,JRR)
END DO
!
@@ -520,7 +521,7 @@ END IF
!
!* 2.3 dissipative heating coeff a t
!
-ZCOEF_DISS(:,:,:) = 1/(ZCP(:,:,:) * ZEXN(:,:,:))
+ZCOEF_DISS(:,:,:) = 1/(ZCP(:,:,:) * ZEXN(:,:,:))
!
!
ZFRAC_ICE(:,:,:) = 0.0
@@ -535,16 +536,26 @@ IF (KRRL >=1) THEN
!
!* 2.5 Lv/Cph/Exn
!
- IF ( KRRI >= 1 ) THEN
+ IF ( KRRI >= 1 ) THEN
ALLOCATE(ZLVOCPEXNM(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)))
ALLOCATE(ZLSOCPEXNM(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)))
ALLOCATE(ZAMOIST_ICE(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)))
ALLOCATE(ZATHETA_ICE(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)))
!
- CALL COMPUTE_FUNCTION_THERMO(XALPW,XBETAW,XGAMW,XLVTT,XCL,ZT,ZEXN,ZCP, &
+!
+ !wc call new functions depending on statnew
+ IF (LSTATNW) THEN
+ CALL COMPUTE_FUNCTION_THERMO_NEW_STAT(XALPW,XBETAW,XGAMW,XLVTT,XCL,ZT,ZEXN,ZCP, &
ZLVOCPEXNM,ZAMOIST,ZATHETA)
- CALL COMPUTE_FUNCTION_THERMO(XALPI,XBETAI,XGAMI,XLSTT,XCI,ZT,ZEXN,ZCP, &
+ CALL COMPUTE_FUNCTION_THERMO_NEW_STAT(XALPI,XBETAI,XGAMI,XLSTT,XCI,ZT,ZEXN,ZCP, &
ZLSOCPEXNM,ZAMOIST_ICE,ZATHETA_ICE)
+ ELSE
+
+ CALL COMPUTE_FUNCTION_THERMO(XALPW,XBETAW,XGAMW,XLVTT,XCL,ZT,ZEXN,ZCP, &
+ ZLVOCPEXNM,ZAMOIST,ZATHETA)
+ CALL COMPUTE_FUNCTION_THERMO(XALPI,XBETAI,XGAMI,XLSTT,XCI,ZT,ZEXN,ZCP, &
+ ZLSOCPEXNM,ZAMOIST_ICE,ZATHETA_ICE)
+ ENDIF
!
WHERE(PRT(:,:,:,2)+PRT(:,:,:,4)>0.0)
ZFRAC_ICE(:,:,:) = PRT(:,:,:,4) / ( PRT(:,:,:,2)+PRT(:,:,:,4) )
@@ -560,8 +571,14 @@ IF (KRRL >=1) THEN
DEALLOCATE(ZAMOIST_ICE)
DEALLOCATE(ZATHETA_ICE)
ELSE
- CALL COMPUTE_FUNCTION_THERMO(XALPW,XBETAW,XGAMW,XLVTT,XCL,ZT,ZEXN,ZCP, &
+ !wc call new stat functions or not
+ IF (LSTATNW) THEN
+ CALL COMPUTE_FUNCTION_THERMO_NEW_STAT(XALPW,XBETAW,XGAMW,XLVTT,XCL,ZT,ZEXN,ZCP, &
ZLOCPEXNM,ZAMOIST,ZATHETA)
+ ELSE
+ CALL COMPUTE_FUNCTION_THERMO(XALPW,XBETAW,XGAMW,XLVTT,XCL,ZT,ZEXN,ZCP, &
+ ZLOCPEXNM,ZAMOIST,ZATHETA)
+ ENDIF
END IF
!
!
@@ -577,7 +594,7 @@ IF (KRRL >=1) THEN
TZFIELD%NDIMS = 3
TZFIELD%LTIMEDEP = .TRUE.
CALL IO_FIELD_WRITE(TPFILE,TZFIELD,ZATHETA)
-!
+!
TZFIELD%CMNHNAME = 'AMOIST'
TZFIELD%CSTDNAME = ''
TZFIELD%CLONGNAME = 'AMOIST'
@@ -609,7 +626,7 @@ IF ( KRRL >= 1 ) THEN
- ZLSOCPEXNM(:,:,:) * PRRS(:,:,:,4)
ELSE
! Rnp at t
- PRT(:,:,:,1) = PRT(:,:,:,1) + PRT(:,:,:,2)
+ PRT(:,:,:,1) = PRT(:,:,:,1) + PRT(:,:,:,2)
PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,2)
! Theta_l at t
PTHLT(:,:,:) = PTHLT(:,:,:) - ZLOCPEXNM(:,:,:) * PRT(:,:,:,2)
@@ -651,7 +668,7 @@ SELECT CASE (HTURBLEN)
ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM,OOCEAN)
!
-!* 3.3 Grey-zone combined RM17 & Deardorff mixing lengths
+!* 3.3 Grey-zone combined RM17 & Deardorff mixing lengths
! --------------------------------------------------
CASE ('ADAP')
@@ -663,7 +680,7 @@ SELECT CASE (HTURBLEN)
CALL DELT(ZLMW,ODZ=.FALSE.)
! The minimum mixing length is chosen between Horizontal grid mesh (not taking into account the vertical grid mesh) and RM17.
! For large horizontal grid meshes, this is equal to RM17
- ! For LES grid meshes, this is equivalent to Deardorff : the base mixing lentgh is the horizontal grid mesh,
+ ! For LES grid meshes, this is equivalent to Deardorff : the base mixing lentgh is the horizontal grid mesh,
! and it is limited by a stability-based length (RM17), as was done in Deardorff length (but taking into account shear as well)
! For grid meshes in the grey zone, then this is the smaller of the two.
ZLM = MIN(ZLM,XCADAP*ZLMW)
@@ -743,7 +760,7 @@ IF (HTURBDIM=="3DIM") THEN
CALL UPDATE_LM(HLBCX,HLBCY,ZLM,ZLEPS)
END IF
!
-!* 3.9 Mixing length correction if immersed walls
+!* 3.9 Mixing length correction if immersed walls
! ------------------------------------------
!
IF (LIBM) THEN
@@ -780,7 +797,7 @@ ZCDUEFF(:,:) =-SQRT ( (PSFU(:,:)**2 + PSFV(:,:)**2) / &
(1.E-60 + ZUSLOPE(:,:)**2 + ZVSLOPE(:,:)**2 ) )
#else
(XMNH_TINY + ZUSLOPE(:,:)**2 + ZVSLOPE(:,:)**2 ) )
-#endif
+#endif
!
!* 4.6 compute the surface tangential fluxes
!
@@ -925,8 +942,8 @@ END IF
!
!Les budgets des termes horizontaux de la turb sont présents dans AROME
! alors que ces termes ne sont pas calculés
-#ifdef REPRO48
-#else
+#ifdef REPRO48
+#else
IF( HTURBDIM == '3DIM' ) THEN
#endif
IF( LBUDGET_U ) CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_U ), 'HTURB', PRUS (:, :, :) )
@@ -1020,11 +1037,11 @@ END IF
#endif
!----------------------------------------------------------------------------
!
-!* 6. EVOLUTION OF THE TKE AND ITS DISSIPATION
+!* 6. EVOLUTION OF THE TKE AND ITS DISSIPATION
! ----------------------------------------
!
-! 6.1 Contribution of mass-flux in the TKE buoyancy production if
-! cloud computation is not statistical
+! 6.1 Contribution of mass-flux in the TKE buoyancy production if
+! cloud computation is not statistical
PTP = PTP + XG / PTHVREF * MZF(PFLXZTHVMF,KKA, KKU, KKL)
IF(PRESENT(PTPMF)) PTPMF=XG / PTHVREF * MZF(PFLXZTHVMF, KKA, KKU, KKL)
@@ -1077,9 +1094,9 @@ ENDIF
! ---------------------------------------------------------
!
IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
-!
+!
! stores the mixing length
-!
+!
TZFIELD%CMNHNAME = 'LM'
TZFIELD%CSTDNAME = ''
TZFIELD%CLONGNAME = 'LM'
@@ -1129,7 +1146,7 @@ IF(PRESENT(PDRUS_TURB)) THEN
PDRUS_TURB = PRUS - PDRUS_TURB
PDRVS_TURB = PRVS - PDRVS_TURB
PDRTHLS_TURB = PRTHLS - PDRTHLS_TURB
- PDRRTS_TURB = PRRS(:,:,:,1) - PDRRTS_TURB
+ PDRRTS_TURB = PRRS(:,:,:,1) - PDRRTS_TURB
PDRSVS_TURB = PRSVS - PDRSVS_TURB
END IF
!----------------------------------------------------------------------------
@@ -1149,7 +1166,7 @@ IF ( KRRL >= 1 ) THEN
DEALLOCATE(ZLVOCPEXNM)
DEALLOCATE(ZLSOCPEXNM)
ELSE
- PRT(:,:,:,1) = PRT(:,:,:,1) - PRT(:,:,:,2)
+ PRT(:,:,:,1) = PRT(:,:,:,1) - PRT(:,:,:,2)
PRRS(:,:,:,1) = PRRS(:,:,:,1) - PRRS(:,:,:,2)
PTHLT(:,:,:) = PTHLT(:,:,:) + ZLOCPEXNM(:,:,:) * PRT(:,:,:,2)
PRTHLS(:,:,:) = PRTHLS(:,:,:) + ZLOCPEXNM(:,:,:) * PRRS(:,:,:,2)
@@ -1333,14 +1350,14 @@ USE MODD_CST
!
IMPLICIT NONE
!
-!* 0.1 Declarations of dummy arguments
+!* 0.1 Declarations of dummy arguments
!
REAL, INTENT(IN) :: PALP,PBETA,PGAM,PLTT,PC
REAL, DIMENSION(:,:,:), INTENT(IN) :: PT,PEXN,PCP
!
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PLOCPEXN
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PAMOIST,PATHETA
-!
+!
!* 0.2 Declarations of local variables
!
REAL :: ZEPS ! XMV / XMD
@@ -1394,6 +1411,84 @@ REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZDRVSATDT
!
IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO',1,ZHOOK_HANDLE)
END SUBROUTINE COMPUTE_FUNCTION_THERMO
+
+! ########################################################################
+ SUBROUTINE COMPUTE_FUNCTION_THERMO_NEW_STAT(PALP,PBETA,PGAM,PLTT,PC,PT,PEXN,PCP,&
+ PLOCPEXN,PAMOIST,PATHETA )
+! ########################################################################
+!!
+!!**** *COMPUTE_FUNCTION_THERMO* routine to compute several thermo functions
+!
+!! AUTHOR
+!! ------
+!!
+!! JP Pinty *LA*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 24/02/03
+!! Modified: Wim de Rooy 06-02-2019
+!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+USE MODD_CST
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments
+!
+REAL, INTENT(IN) :: PALP,PBETA,PGAM,PLTT,PC
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PT,PEXN,PCP
+!
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PLOCPEXN
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PAMOIST,PATHETA
+!
+!* 0.2 Declarations of local variables
+!
+REAL :: ZEPS ! XMV / XMD
+REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZRVSAT
+REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZDRVSATDT
+!
+!-------------------------------------------------------------------------------
+!
+ REAL(KIND=JPRB) :: ZHOOK_HANDLE
+ IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO_NEW_STAT',0,ZHOOK_HANDLE)
+ ZEPS = XMV / XMD
+!
+!* 1.1 Lv/Cph at t
+!
+ PLOCPEXN(:,:,:) = ( PLTT + (XCPV-PC) * (PT(:,:,:)-XTT) ) / PCP(:,:,:)
+!
+!* 1.2 Saturation vapor pressure at t
+!
+ ZRVSAT(:,:,:) = EXP( PALP - PBETA/PT(:,:,:) - PGAM*ALOG( PT(:,:,:) ) )
+!
+!* 1.3 saturation mixing ratio at t
+!
+ ZRVSAT(:,:,:) = ZRVSAT(:,:,:) * ZEPS / ( PPABST(:,:,:) - ZRVSAT(:,:,:) )
+!
+!* 1.4 compute the saturation mixing ratio derivative (rvs')
+!
+ ZDRVSATDT(:,:,:) = ( PBETA / PT(:,:,:) - PGAM ) / PT(:,:,:) &
+ * ZRVSAT(:,:,:) * ( 1. + ZRVSAT(:,:,:) / ZEPS )
+!
+!* 1.5 compute Amoist
+!
+ PAMOIST(:,:,:)= 1.0 / ( 1.0 + ZDRVSATDT(:,:,:) * PLOCPEXN(:,:,:) )
+!
+!* 1.6 compute Atheta
+!
+ PATHETA(:,:,:)= PAMOIST(:,:,:) * PEXN(:,:,:) * ZDRVSATDT(:,:,:)
+!
+!* 1.7 Lv/Cph/Exner at t-1
+!
+ PLOCPEXN(:,:,:) = PLOCPEXN(:,:,:) / PEXN(:,:,:)
+!
+IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO_NEW_STAT',1,ZHOOK_HANDLE)
+END SUBROUTINE COMPUTE_FUNCTION_THERMO_NEW_STAT
+
!
! ####################
SUBROUTINE DELT(PLM,ODZ)
@@ -1415,7 +1510,7 @@ END SUBROUTINE COMPUTE_FUNCTION_THERMO
!* 0. DECLARATIONS
! ------------
!
-!* 0.1 Declarations of dummy arguments
+!* 0.1 Declarations of dummy arguments
!
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PLM
LOGICAL, INTENT(IN) :: ODZ
@@ -1437,7 +1532,7 @@ IF (ODZ) THEN
PLM(:,:,KKA) = PZZ(:,:,IKB) - PZZ(:,:,KKA)
IF ( HTURBDIM /= '1DIM' ) THEN ! 3D turbulence scheme
IF ( L2D) THEN
- PLM(:,:,:) = SQRT( PLM(:,:,:)*MXF(PDXX(:,:,:)) )
+ PLM(:,:,:) = SQRT( PLM(:,:,:)*MXF(PDXX(:,:,:)) )
ELSE
PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.)
END IF
@@ -1454,7 +1549,7 @@ ELSE
END IF
END IF
!
-! mixing length limited by the distance normal to the surface
+! mixing length limited by the distance normal to the surface
! (with the same factor as for BL89)
!
IF (.NOT. ORMC01) THEN
@@ -1481,7 +1576,7 @@ IF (.NOT. ORMC01) THEN
EXIT
ENDIF
END DO
- ENDIF
+ ENDIF
END DO
END DO
END IF
@@ -1514,7 +1609,7 @@ END SUBROUTINE DELT
!* 0. DECLARATIONS
! ------------
!
-!* 0.1 Declarations of dummy arguments
+!* 0.1 Declarations of dummy arguments
!
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PLM
!
@@ -1526,7 +1621,7 @@ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZWORK2D
!
REAL, DIMENSION(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)) :: &
ZDTHLDZ,ZDRTDZ, &!dtheta_l/dz, drt_dz used for computing the stablity
-! ! criterion
+! ! criterion
ZETHETA,ZEMOIST !coef ETHETA and EMOIST
!----------------------------------------------------------------------------
!
@@ -1593,7 +1688,7 @@ ELSE! For dry atmos or unsalted ocean runs
END DO
END DO
END IF
-! special case near the surface
+! special case near the surface
ZDTHLDZ(:,:,IKB)=(PTHLT(:,:,IKB+KKL)-PTHLT(:,:,IKB))/PDZZ(:,:,IKB+KKL)
! For dry simulations
IF (KRR>0) THEN
@@ -1639,7 +1734,7 @@ IF (.NOT. ORMC01) THEN
EXIT
ENDIF
END DO
- ENDIF
+ ENDIF
END DO
END DO
END IF
@@ -1707,7 +1802,7 @@ REAL :: ZCOEF_AMPL_CEI_NUL! Ordonnate at the origin of the
! amplification straight line
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZCOEF_AMPL
! Amplification coefficient of the mixing length
- ! when the instability criterium is verified
+ ! when the instability criterium is verified
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZLM_CLOUD
! Turbulent mixing length in the clouds
!
@@ -1718,7 +1813,7 @@ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZLM_CLOUD
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB:CLOUD_MODIF_LM',0,ZHOOK_HANDLE)
-ZPENTE = ( PCOEF_AMPL_SAT - 1. ) / ( PCEI_MAX - PCEI_MIN )
+ZPENTE = ( PCOEF_AMPL_SAT - 1. ) / ( PCEI_MAX - PCEI_MIN )
ZCOEF_AMPL_CEI_NUL = 1. - ZPENTE * PCEI_MIN
!
ZCOEF_AMPL(:,:,:) = 1.
@@ -1735,7 +1830,7 @@ WHERE ( PCEI(:,:,:)>=PCEI_MAX ) ZCOEF_AMPL(:,:,:)=PCOEF_AMPL_SAT
!
WHERE ( PCEI(:,:,:) < PCEI_MAX .AND. &
PCEI(:,:,:) > PCEI_MIN ) &
- ZCOEF_AMPL(:,:,:) = ZPENTE * PCEI(:,:,:) + ZCOEF_AMPL_CEI_NUL
+ ZCOEF_AMPL(:,:,:) = ZPENTE * PCEI(:,:,:) + ZCOEF_AMPL_CEI_NUL
!
!
!* 3. CALCULATION OF THE MIXING LENGTH IN CLOUDS
@@ -1824,4 +1919,4 @@ ENDIF
IF (LHOOK) CALL DR_HOOK('TURB:CLOUD_MODIF_LM',1,ZHOOK_HANDLE)
END SUBROUTINE CLOUD_MODIF_LM
!
-END SUBROUTINE TURB
+END SUBROUTINE TURB