Skip to content
Snippets Groups Projects
Commit 63323d78 authored by RODIER Quentin's avatar RODIER Quentin
Browse files

Quentin 08/01/2019: correction of the TKE dissipation constant with RM17...

Quentin 08/01/2019: correction of the TKE dissipation constant with RM17 mixing length and add RM17 in shallow convection scheme
- ini_cturb.f90 + read_exsegn.f90 : XCED constant is replaced by XCEDIS which depends on BL89 or RM17 mixing length (see Rodier et al. 2017)
- shallow_mf.f90 + all compute_*.f90 : support of RM17 mixing length in shallow convection
- bl89.f90 : add comment
parent 170b3f07
No related branches found
No related tags found
No related merge requests found
......@@ -71,6 +71,7 @@ END MODULE MODI_BL89
!! 2012-02 (Y. Seity) add possibility to run with
!! reversed vertical levels
!! Philippe 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q
!! 01/2019 (Q. Rodier) support for RM17 mixing length
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
......@@ -223,12 +224,12 @@ ZHLVPT(:,KKA) = ZVPT(:,KKA)
!
!* 3. loop on model levels
! --------------------
!
DO JK=IKTB,IKTE
!
!-------------------------------------------------------------------------------
!
!
!* 4. mixing length for a downwards displacement
! ------------------------------------------
ZINTE(:)=ZTKEM(:,JK)
......
......@@ -10,7 +10,7 @@ INTERFACE
! ###################################################################
SUBROUTINE COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ2D, &
PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,PLWORK)
PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,PSHEAR,PLWORK)
! ###################################################################
!* 1.1 Declaration of Arguments
......@@ -29,6 +29,7 @@ LOGICAL, INTENT(IN) :: OUPORDN ! switch to compute upward
! downward (false) mixing length
LOGICAL, INTENT(IN) :: OFLUX ! Computation must be done from flux level
REAL, DIMENSION(:), INTENT(OUT) :: PLWORK ! Resulting mixing length
REAL, DIMENSION(:,:), INTENT(IN) :: PSHEAR ! vertical wind shear for RM17 mixing length
END SUBROUTINE COMPUTE_BL89_ML
......@@ -37,7 +38,7 @@ END INTERFACE
END MODULE MODI_COMPUTE_BL89_ML
! ######spl
SUBROUTINE COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ2D, &
PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,PLWORK)
PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,PSHEAR,PLWORK)
! ###################################################################
!!
......@@ -53,6 +54,7 @@ END MODULE MODI_COMPUTE_BL89_ML
!! Original 19/01/06
!! S. Riette Jan 2012: support for both order of vertical levels and cleaning
!! R.Honnert Oct 2016 : Update with AROME
!! Q.Rodier 01/2019 : support RM17 mixing length as in bl89.f90
!!
!-------------------------------------------------------------------------------
!
......@@ -92,6 +94,7 @@ LOGICAL, INTENT(IN) :: OUPORDN ! switch to compute upward
! downward (false) mixing length
LOGICAL, INTENT(IN) :: OFLUX ! Computation must be done from flux level
REAL, DIMENSION(:), INTENT(OUT) :: PLWORK ! Resulting mixing length
REAL, DIMENSION(:,:), INTENT(IN) :: PSHEAR ! vertical wind shear for RM17 mixing length
! 0.2 Local variable
!
......@@ -144,21 +147,24 @@ IF (OUPORDN.EQV..TRUE.) THEN
ZTEST0=0.5+SIGN(0.5,ZINTE(J1D)) ! test if there's energy to consume
! Energy consumed if parcel cross the entire layer
ZPOTE(J1D) = ZTEST0*(PG_O_THVREF(J1D) * &
(0.5*(ZHLVPT(J1D,KK)+ PVPT(J1D,KK)) - ZVPT_DEP(J1D))) * &
(0.5*(ZHLVPT(J1D,KK)+ PVPT(J1D,KK)) - ZVPT_DEP(J1D)) + &
XRM17*PSHEAR(J1D,KK)*SQRT(ABS(PTKEM_DEP(J1D)))) * &
PDZZ2D(J1D,KK)*0.5
! Test if it rests some energy to consume
ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE(J1D))
! Length travelled by parcel if it rests energy to consume
ZLWORK1(J1D)=PDZZ2D(J1D,KK)*0.5
! Lenght travelled by parcel to nullify energy
ZLWORK2(J1D)= ( - PG_O_THVREF(J1D) * &
ZLWORK2(J1D)= ( - PG_O_THVREF(J1D) * &
( ZHLVPT(J1D,KK) - ZVPT_DEP(J1D) ) &
- XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) &
+ SQRT (ABS( &
( PG_O_THVREF(J1D) * (ZHLVPT(J1D,KK) - ZVPT_DEP(J1D)) )**2 &
(XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) + &
PG_O_THVREF(J1D) * (ZHLVPT(J1D,KK) - ZVPT_DEP(J1D)) )**2 &
+ 2. * ZINTE(J1D) * PG_O_THVREF(J1D) &
* ZDELTVPT(J1D,KK) / PDZZ2D(J1D,KK) )) ) / &
( PG_O_THVREF(J1D) * ZDELTVPT(J1D,KK) / PDZZ2D(J1D,KK) )
! Effective length travelled by parcel
( PG_O_THVREF(J1D) * ZDELTVPT(J1D,KK) / PDZZ2D(J1D,KK) )
! Effective length travelled by parcel
PLWORK(J1D)=PLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1(J1D)+ &
(1-ZTEST)*ZLWORK2(J1D))
! Rest of energy to consume
......@@ -174,19 +180,22 @@ IF (OUPORDN.EQV..TRUE.) THEN
DO J1D=1,IIJU
ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
ZPOTE(J1D) = ZTEST0*(PG_O_THVREF(J1D) * &
(ZHLVPT(J1D,JKK) - ZVPT_DEP(J1D))) * PDZZ2D(J1D,JKK) !particle keeps its temperature
(ZHLVPT(J1D,JKK) - ZVPT_DEP(J1D)) &
+ XRM17*PSHEAR(J1D,JKK)*SQRT(ABS(PTKEM_DEP(J1D))))* PDZZ2D(J1D,JKK)
ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE(J1D))
ZTESTM=ZTESTM+ZTEST0
ZLWORK1(J1D)=PDZZ2D(J1D,JKK)
!ZLWORK2 jump of the last reached level
ZLWORK2(J1D)= ( - PG_O_THVREF(J1D) * &
( PVPT(J1D,JKK-KKL) - ZVPT_DEP(J1D) ) &
+ SQRT (ABS( &
( PG_O_THVREF(J1D) * (PVPT(J1D,JKK-KKL) - ZVPT_DEP(J1D)) )**2 &
- XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) &
+ SQRT (ABS( &
(XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) + &
PG_O_THVREF(J1D) * (PVPT(J1D,JKK-KKL) - ZVPT_DEP(J1D)) )**2 &
+ 2. * ZINTE(J1D) * PG_O_THVREF(J1D) &
* ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )) ) / &
( PG_O_THVREF(J1D) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
!
!
PLWORK(J1D)=PLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1(J1D)+ &
(1-ZTEST)*ZLWORK2(J1D))
ZINTE(J1D) = ZINTE(J1D) - ZPOTE(J1D)
......@@ -214,19 +223,22 @@ IF (OUPORDN.EQV..FALSE.) THEN
ZTESTM=0
DO J1D=1,IIJU
ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
ZPOTE(J1D) = -ZTEST0*(PG_O_THVREF(J1D) * &
(ZHLVPT(J1D,JKK) - PVPT(J1D,KK))) * PDZZ2D(J1D,JKK) !particle keeps its temperature
ZPOTE(J1D) = ZTEST0*(-PG_O_THVREF(J1D) * &
(ZHLVPT(J1D,JKK) - PVPT(J1D,KK)) &
+ XRM17*PSHEAR(J1D,JKK)*SQRT(ABS(PTKEM_DEP(J1D))))* PDZZ2D(J1D,JKK)
ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE(J1D))
ZTESTM=ZTESTM+ZTEST0
ZLWORK1(J1D)=PDZZ2D(J1D,JKK)
ZLWORK2(J1D)= ( + PG_O_THVREF(J1D) * &
ZLWORK2(J1D)= ( + PG_O_THVREF(J1D) * &
( PVPT(J1D,JKK) - PVPT(J1D,KK) ) &
-XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) &
+ SQRT (ABS( &
( PG_O_THVREF(J1D) * (PVPT(J1D,JKK) - PVPT(J1D,KK)) )**2 &
(XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) - &
PG_O_THVREF(J1D) * (PVPT(J1D,JKK) - PVPT(J1D,KK)) )**2 &
+ 2. * ZINTE(J1D) * PG_O_THVREF(J1D) &
* ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )) ) / &
( PG_O_THVREF(J1D) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
!
!
PLWORK(J1D)=PLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1(J1D)+ &
(1-ZTEST)*ZLWORK2(J1D))
ZINTE(J1D) = ZINTE(J1D) - ZPOTE(J1D)
......@@ -234,5 +246,5 @@ IF (OUPORDN.EQV..FALSE.) THEN
ENDIF
END DO
ENDIF
END SUBROUTINE COMPUTE_BL89_ML
......@@ -134,6 +134,7 @@ END MODULE MODI_COMPUTE_UPDRAFT
!! V.Masson, C.Lac : 02/2011 : SV_UP initialized by a non-zero value
!! S. Riette Apr 2013: improvement of continuity at the condensation level
!! R.Honnert Oct 2016 : Add ZSURF and Update with AROME
!! Q.Rodier 01/2019 : support RM17 mixing length
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
......@@ -141,6 +142,7 @@ END MODULE MODI_COMPUTE_UPDRAFT
!
USE MODD_CST
USE MODD_PARAM_MFSHALL_n
USE MODD_TURB_n, ONLY : CTURBLEN
USE MODI_COMPUTE_ENTR_DETR
USE MODI_TH_R_FROM_THL_RT_1D
......@@ -258,7 +260,7 @@ REAL :: ZDEPTH_MAX1, ZDEPTH_MAX2 ! control auto-extinction process
REAL :: ZTMAX,ZRMAX ! control value
REAL, DIMENSION(SIZE(PTHM,1)) :: ZSURF
REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZSHEAR,ZDUDZ,ZDVDZ ! vertical wind shear
! Thresholds for the perturbation of
! theta_l and r_t at the first level of the updraft
ZTMAX=2.0
......@@ -377,8 +379,16 @@ IF (OENTR_DETR) THEN
! compute L_up
GLMIX=.TRUE.
ZTKEM_F(:,KKB)=0.
CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB),ZTHVM,KKB,GLMIX,.FALSE.,ZLUP)
!
IF(CTURBLEN=='RM17') THEN
ZDUDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PUM,PDZZ))
ZDVDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PVM,PDZZ))
ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
ELSE
ZSHEAR = 0. !no shear in bl89 mixing length
END IF
!
CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB),ZTHVM,KKB,GLMIX,.FALSE.,ZSHEAR,ZLUP)
ZLUP(:)=MAX(ZLUP(:),1.E-10)
! Compute Buoyancy flux at the ground
......
......@@ -134,6 +134,7 @@ END MODULE MODI_COMPUTE_UPDRAFT_HRIO! ######spl
!! S. Riette may 2011: ice added, interface modified
!! S. Riette Jan 2012: support for both order of vertical levels
!! V.Masson, C.Lac : 02/2011 : SV_UP initialized by a non-zero value
!! Q.Rodier 01/2019 : support RM17 mixing length
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
......@@ -147,6 +148,7 @@ USE MODD_PARAM_MFSHALL_n, ONLY : XPRES_UV,XALP_PERT,XCMF,XFRAC_UP_MAX,XA1,XB,&
! XC,XBETA1
USE MODD_GRID_n, ONLY : XDXHAT, XDYHAT
USE MODD_BLANK
USE MODD_TURB_n, ONLY :CTURBLEN
!USE MODI_COMPUTE_ENTR_DETR
USE MODI_TH_R_FROM_THL_RT_1D
......@@ -275,6 +277,7 @@ REAL :: XFRAC_LIM ! surface maximale du thermique
REAL :: ZTMAX,ZRMAX, ZEPS ! control value
REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZSHEAR,ZDUDZ,ZDVDZ ! vertical wind shear
! pour le calcul de la resolution normalise
REAL, DIMENSION(SIZE(PTHM,1)) :: ZA1, ZRESOL_NORM
REAL :: ZRESOL_GRID
......@@ -446,9 +449,15 @@ IF (OENTR_DETR) THEN
! compute L_up
GLMIX=.TRUE.
ZTKEM_F(:,KKB)=0.
IF(CTURBLEN=='RM17') THEN
ZDUDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PUM,PDZZ))
ZDVDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PVM,PDZZ))
ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
ELSE
ZSHEAR = 0. !no shear in bl89 mixing length
END IF
CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB), &
ZTHVM_F,KKB,GLMIX,.TRUE.,ZLUP)
ZTHVM_F,KKB,GLMIX,.TRUE.,ZSHEAR,ZLUP)
ZLUP(:)=MAX(ZLUP(:),1.E-10)
! Compute Buoyancy flux at the ground
......
......@@ -124,6 +124,7 @@ SUBROUTINE COMPUTE_UPDRAFT_RHCJ10(KKA,KKB,KKE,KKU,KKL,HFRAC_ICE, &
!! ------
!! Y. Bouteloup (2012)
!! R. Honert Janv 2013 ==> corection of some bugs
!! Q.Rodier 01/2019 : support RM17 mixing length
!! --------------------------------------------------------------------------
! WARNING ==> This updraft is not yet ready to use scalar variables
......@@ -133,7 +134,7 @@ SUBROUTINE COMPUTE_UPDRAFT_RHCJ10(KKA,KKB,KKE,KKU,KKL,HFRAC_ICE, &
USE MODD_CST
USE MODD_PARAM_MFSHALL_n
USE MODD_TURB_n, ONLY : CTURBLEN
USE MODI_COMPUTE_ENTR_DETR
USE MODI_TH_R_FROM_THL_RT_1D
USE MODI_SHUMAN_MF
......@@ -252,6 +253,7 @@ REAL :: ZDEPTH_MAX1, ZDEPTH_MAX2 ! control auto-extinction process
REAL :: ZTMAX,ZRMAX, ZEPS ! control value
REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZSHEAR,ZDUDZ,ZDVDZ ! vertical wind shear
! Thresholds for the perturbation of
! theta_l and r_t at the first level of the updraft
......@@ -389,9 +391,17 @@ ENDDO
! compute L_up
GLMIX=.TRUE.
ZTKEM_F(:,KKB)=0.
!
IF(CTURBLEN=='RM17') THEN
ZDUDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PUM,PDZZ))
ZDVDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PVM,PDZZ))
ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
ELSE
ZSHEAR = 0. !no shear in bl89 mixing length
END IF
!
CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB), &
ZTHVM_F,KKB,GLMIX,.TRUE.,ZLUP)
ZTHVM_F,KKB,GLMIX,.TRUE.,ZSHEAR,ZLUP)
ZLUP(:)=MAX(ZLUP(:),1.E-10)
! Compute Buoyancy flux at the ground
......
......@@ -63,6 +63,7 @@ END MODULE MODI_INI_CTURB
!! P.Jabouille 20/10/99 XCET=0.4
!! V.Masson 13/11/02 XALPSBL and XASBL
!! 05/06 Remove KEPS
!! Q.Rodier 01/19 XCED replaced by XCEDIS in read_exsegn.f90 and ini_modeln.f90
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
......@@ -80,8 +81,10 @@ IMPLICIT NONE
!
! 1.1 Constant for dissipation of Tke
!
!XCED is now replaced by XCEDIS
!XCED = 0.70
XCED = 0.84
!XCED = 0.84
!
! Redelsperger-Sommeria (1981) = 0.70
! Schmidt-Schumann (1989) = 0.845
! Cheng-Canuto-Howard (2002) = 0.845
......
......@@ -288,6 +288,7 @@ END MODULE MODI_READ_EXSEG_n
!! Q.Libois 02/2018 ECRAD
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
!! Modification 07/2017 (V. Vionnet) add blowing snow scheme
!! Modification 01/2019 (Q. Rodier) define XCEDIS depending on BL89 or RM17 mixing length
!!------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
......@@ -1414,7 +1415,12 @@ ELSE
CGETCLDFR = 'SKIP'
END IF
END IF
!
IF(CTURBLEN=='RM17') THEN
XCEDIS=0.34
ELSE
XCEDIS=0.84
END IF
!
!* 3.3 Moist turbulence
!
......
......@@ -165,6 +165,7 @@ END MODULE MODI_SHALLOW_MF
!! R.Honnert 10/2016 : SURF=gray zone initilisation + EDKF
!! R.Honnert 10/2016 : Update with Arome
!! Philippe Wautelet 28/05/2018: corrected truncated integer division (2/3 -> 2./3.)
!! Q.Rodier 01/2019 : support RM17 mixing length
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
......@@ -173,6 +174,7 @@ END MODULE MODI_SHALLOW_MF
USE MODD_CST
USE MODD_PARAMETERS, ONLY: JPVEXT
USE MODD_PARAM_MFSHALL_n
USE MODD_TURB_n, ONLY: CTURBLEN
USE MODI_THL_RT_FROM_TH_R_MF
USE MODI_COMPUTE_UPDRAFT
......@@ -184,7 +186,7 @@ USE MODI_MF_TURB_EXPL
USE MODI_MF_TURB_GREYZONE
USE MODI_COMPUTE_MF_CLOUD
USE MODI_COMPUTE_FRAC_ICE
!
USE MODI_SHUMAN_MF
!
USE MODI_COMPUTE_BL89_ML
USE MODD_GRID_n, ONLY : XDXHAT, XDYHAT
......@@ -282,6 +284,7 @@ REAL, DIMENSION(SIZE(PSVM,1),SIZE(PSVM,2),SIZE(PSVM,3)) :: &
REAL, DIMENSION(SIZE(PTHM,1)) :: ZDEPTH ! Deepness of cloud
REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZFRAC_ICE_UP ! liquid/solid fraction in updraft
REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZRSAT_UP ! Rsat in updraft
REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZSHEAR,ZDUDZ,ZDVDZ !vertical wind shear
LOGICAL :: GENTR_DETR ! flag to recompute entrainment, detrainment and mass flux
INTEGER :: IKB ! near ground physical index
......@@ -460,7 +463,14 @@ ENDIF
ENDDO
ZG_O_THVREF=XG/PTHVREF
GLMIX=.TRUE.
CALL COMPUTE_BL89_ML(KKA,IKB,IKE,KKU,KKL,PDZZ,PTKEM(:,IKB) ,ZG_O_THVREF(:,IKB),ZTHVM,IKB,GLMIX,.TRUE.,ZLUP)
IF(CTURBLEN=='RM17') THEN
ZDUDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PUM,PDZZ))
ZDVDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PVM,PDZZ))
ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ)
ELSE
ZSHEAR = 0. !no shear in bl89 mixing length
END IF
CALL COMPUTE_BL89_ML(KKA,IKB,IKE,KKU,KKL,PDZZ,PTKEM(:,IKB) ,ZG_O_THVREF(:,IKB),ZTHVM,IKB,GLMIX,.TRUE.,ZSHEAR,ZLUP)
!! calcul de Dx/(h+hc)
DO JI=1,SIZE(XDXHAT)
DO JJ=1,SIZE(XDYHAT)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment