diff --git a/src/MNH/bl89.f90 b/src/MNH/bl89.f90
index bb81821ef0dbb6c51471c03d40d0d55d6c4477bc..6835f449241cc58aed7328184012f2d70bd5036c 100644
--- a/src/MNH/bl89.f90
+++ b/src/MNH/bl89.f90
@@ -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)
diff --git a/src/MNH/compute_bl89_ml.f90 b/src/MNH/compute_bl89_ml.f90
index 4c5124602b87fefe4ceddf7175abe469eeef5647..ac8986c3df7b608442bd86a414a8496e3c99e9f1 100644
--- a/src/MNH/compute_bl89_ml.f90
+++ b/src/MNH/compute_bl89_ml.f90
@@ -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
diff --git a/src/MNH/compute_updraft.f90 b/src/MNH/compute_updraft.f90
index e4165877555f575709ee19b94dc96dd1e43ee46a..ba78d05a1fdaec9883a5eb1685fcdd3935d4c3cd 100644
--- a/src/MNH/compute_updraft.f90
+++ b/src/MNH/compute_updraft.f90
@@ -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
diff --git a/src/MNH/compute_updraft_hrio.f90 b/src/MNH/compute_updraft_hrio.f90
index 3466441418fdbdd0b983c56abc4b5e807aba9549..6ec162b20d6bd4c14587325a2e5e4c266ca4dc19 100644
--- a/src/MNH/compute_updraft_hrio.f90
+++ b/src/MNH/compute_updraft_hrio.f90
@@ -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 normalisée
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
diff --git a/src/MNH/compute_updraft_rhcj10.f90 b/src/MNH/compute_updraft_rhcj10.f90
index f062dcf0ad9d0d6738b16c9fcf1364bb6d2165e1..554dc54f9d8f56fab536ff264a5a989b1d983f23 100644
--- a/src/MNH/compute_updraft_rhcj10.f90
+++ b/src/MNH/compute_updraft_rhcj10.f90
@@ -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
diff --git a/src/MNH/ini_cturb.f90 b/src/MNH/ini_cturb.f90
index 3e7bc3275572eb591a2459ed12471bffb4824d30..090f8ad30de5e7f751e058004b61a318196f9882 100644
--- a/src/MNH/ini_cturb.f90
+++ b/src/MNH/ini_cturb.f90
@@ -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
diff --git a/src/MNH/prep_pgd.f90 b/src/MNH/prep_pgd.f90
index fefe3a572160f3fbd8d14f4b4791eb4a43ced4c8..d7b5899bc4db7c5bfbe2970a1531c02c9b30e40c 100644
--- a/src/MNH/prep_pgd.f90
+++ b/src/MNH/prep_pgd.f90
@@ -72,6 +72,8 @@
!! 10/2016 (S.Faroux S.Bielli) correction for NHALO=0
!! 01/2018 (G.Delautier) SURFEX 8.1
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Q. Rodier 01/2019 : add a new filtering for very high slopes in NAM_ZSFILTER
+!!
!----------------------------------------------------------------------------
!
!* 0. DECLARATION
@@ -131,8 +133,9 @@ CHARACTER(LEN=28) :: YDAD =' ' ! name of dad of input FM file
CHARACTER(LEN=28) :: CPGDFILE ='PGDFILE' ! name of the output file
CHARACTER(LEN=100) :: YMSG
INTEGER :: NZSFILTER=1 ! number of iteration for filter for fine orography
-LOGICAL :: LHSLOP=.FALSE. ! filtering of slopes higher than XHSLOP
-REAL :: XHSLOP=1.2 ! if LHSLOP filtering of slopes higher than XHSLOP
+INTEGER :: NLOCZSFILTER=3 ! number of iteration for filter of local fine orography
+LOGICAL :: LHSLOP=.FALSE. ! filtering of local slopes higher than XHSLOP
+REAL :: XHSLOP=1.0 ! slopes where the local fine filtering is applied
INTEGER :: NSLEVE =12 ! number of iteration for filter for smooth orography
REAL :: XSMOOTH_ZS = XUNDEF ! optional uniform smooth orography for SLEVE coordinate
REAL, DIMENSION(:,:),ALLOCATABLE :: ZWORK ! work array for lat and lon reshape
@@ -144,7 +147,7 @@ TYPE(TFILEDATA),POINTER :: TZFILE => NULL()
TYPE(TFILEDATA),POINTER :: TZNMLFILE => NULL() ! Namelist file
!
NAMELIST/NAM_PGDFILE/CPGDFILE, NHALO
-NAMELIST/NAM_ZSFILTER/NZSFILTER,LHSLOP,XHSLOP
+NAMELIST/NAM_ZSFILTER/NZSFILTER,NLOCZSFILTER,LHSLOP,XHSLOP
NAMELIST/NAM_SLEVE/NSLEVE, XSMOOTH_ZS
NAMELIST/NAM_CONF_PGD/JPHEXT, NHALO_MNH
!------------------------------------------------------------------------------
@@ -278,7 +281,7 @@ NULLIFY(TFILE_SURFEX) !Probably not necessary
!
!* 4. Computes and writes smooth orography for SLEVE coordinate
! ---------------------------------------------------------
-CALL ZSMT_PGD(TZFILE,NZSFILTER,NSLEVE,XSMOOTH_ZS)
+CALL ZSMT_PGD(TZFILE,NZSFILTER,NSLEVE,NLOCZSFILTER,LHSLOP,XHSLOP,XSMOOTH_ZS)
!
IF (.NOT.LCARTESIAN) THEN
!!!! WRITE LAT and LON
diff --git a/src/MNH/read_exsegn.f90 b/src/MNH/read_exsegn.f90
index 28c1e5457fce7826746d6649105bd9a79d5ed8d0..b3340610750f787598f88977e6bf9102d8d50935 100644
--- a/src/MNH/read_exsegn.f90
+++ b/src/MNH/read_exsegn.f90
@@ -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
!
diff --git a/src/MNH/shallow_mf.f90 b/src/MNH/shallow_mf.f90
index 049954b62490fd8a03da54a5359775ba9bd1aa32..2a9e824398bd70a650f87ba4d41237ee0642ff54 100644
--- a/src/MNH/shallow_mf.f90
+++ b/src/MNH/shallow_mf.f90
@@ -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)
diff --git a/src/MNH/zsmt_pgd.f90 b/src/MNH/zsmt_pgd.f90
index 0de2cd8a28051af6f2cdf8b303df64aced7ddad8..14a4be23bbce0a159b70a471aad86fc578fcb36d 100644
--- a/src/MNH/zsmt_pgd.f90
+++ b/src/MNH/zsmt_pgd.f90
@@ -8,14 +8,17 @@
!
INTERFACE
!
- SUBROUTINE ZSMT_PGD(TPFILE,KZSFILTER,KSLEVE,PSMOOTH_ZS)
+ SUBROUTINE ZSMT_PGD(TPFILE,KZSFILTER,KSLEVE,KLOCZSFILTER,OHSLOP,PHSLOP,PSMOOTH_ZS)
!
USE MODD_IO_ll, ONLY : TFILEDATA
!
-TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! File characteristics
-INTEGER, INTENT(IN) :: KZSFILTER ! number of iterations for fine orography
-INTEGER, INTENT(IN) :: KSLEVE ! number of iterations
-REAL, INTENT(IN) :: PSMOOTH_ZS ! optional uniform smooth orography for SLEVE coordinate
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! File characteristics
+INTEGER, INTENT(IN) :: KZSFILTER ! number of iterations for fine orography
+INTEGER, INTENT(IN) :: KSLEVE ! number of iterations
+INTEGER, INTENT(IN) :: KLOCZSFILTER ! number of iteration for filter of local fine orography
+LOGICAL, INTENT(IN) :: OHSLOP ! filtering of local slopes higher than XHSLOP
+REAL, INTENT(IN) :: PHSLOP ! slopes where the local fine filtering is applied
+REAL, INTENT(IN) :: PSMOOTH_ZS ! optional uniform smooth orography for SLEVE coordinate
!
END SUBROUTINE ZSMT_PGD
!
@@ -26,7 +29,7 @@ END MODULE MODI_ZSMT_PGD
!
!
! #############################
- SUBROUTINE ZSMT_PGD(TPFILE,KZSFILTER,KSLEVE,PSMOOTH_ZS)
+ SUBROUTINE ZSMT_PGD(TPFILE,KZSFILTER,KSLEVE,KLOCZSFILTER,OHSLOP,PHSLOP,PSMOOTH_ZS)
! #############################
!
!!**** *ZSMT_PGD* computes smoothed orography for SLEVE coordinate
@@ -56,7 +59,8 @@ END MODULE MODI_ZSMT_PGD
!! -------------
!! Original nov 2005
!! J.Escobar 23/06/2015 : correction for JPHEXT<>1
-!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Q. Rodier 01/2019 : add a new filtering for very high slopes (applied locally)
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -70,15 +74,20 @@ USE MODE_FMREAD
USE MODE_FMWRIT
USE MODE_ll , ONLY : GET_DIM_EXT_ll , ADD2DFIELD_ll , CLEANLIST_ll , UPDATE_HALO_ll
USE MODD_ARGSLIST_ll, ONLY : LIST_ll
+USE MODI_SUM_ll
+USE MODE_FIELD, ONLY: TFIELDDATA, TYPEREAL
!
IMPLICIT NONE
!
!* 0.1 declarations of arguments
!
-TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! File characteristics
-INTEGER, INTENT(IN) :: KZSFILTER ! number of iterations for fine orography
-INTEGER, INTENT(IN) :: KSLEVE ! number of iterations
-REAL, INTENT(IN) :: PSMOOTH_ZS ! optional uniform smooth orography for SLEVE coordinate
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! File characteristics
+INTEGER, INTENT(IN) :: KZSFILTER ! number of iterations for fine orography
+INTEGER, INTENT(IN) :: KSLEVE ! number of iterations
+INTEGER, INTENT(IN) :: KLOCZSFILTER ! number of iteration for filter of local fine orography
+LOGICAL, INTENT(IN) :: OHSLOP ! filtering of local slopes higher than XHSLOP
+REAL, INTENT(IN) :: PHSLOP ! slopes where the local fine filtering is applied
+REAL, INTENT(IN) :: PSMOOTH_ZS ! optional uniform smooth orography for SLEVE coordinate
!
!
!* 0.2 declarations of local variables
@@ -88,8 +97,6 @@ INTEGER :: JN ! loop counter on iterations
INTEGER :: JI ! loop counter on X coordinate
INTEGER :: JJ ! loop counter on Y coordinate
!
-INTEGER :: IIMAX ! number of physical points in X direction
-INTEGER :: IJMAX ! number of physical points in Y direction
INTEGER :: IIU ! number of points in X direction
INTEGER :: IJU ! number of points in Y direction
!
@@ -98,11 +105,22 @@ REAL, DIMENSION(:,:), ALLOCATABLE :: ZFINE_ZS ! smoothed fine orography
REAL, DIMENSION(:,:), ALLOCATABLE :: ZSLEVE_ZS ! smooth orography for sleve coordinate
REAL, DIMENSION(:,:), ALLOCATABLE :: ZZS ! orography at previous iteration
REAL, DIMENSION(:,:), ALLOCATABLE :: ZMASK ! sea mask
+!
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSLOPEX ! terrain slope along x (flux pt)
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSLOPEY ! terrain slope along y (flux pt)
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSLOPEMX ! terrain slope along x (max. of flux, at mass pt)
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSLOPEMY ! terrain slope along y (max. of flux, at mass pt)
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOFIL ! coefficient for filtering high slopes
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSMOOTH_ZSINI ! initial orography before slope filtering
+!
INTEGER :: JIM, JIP, JJM, JJP
TYPE(LIST_ll) , POINTER :: THALO_ll => NULL() ! halo
INTEGER :: INFO_ll ! error return code
INTEGER :: IIB,IIE,IJB,IJE
+REAL, DIMENSION(:), ALLOCATABLE :: ZXHAT
+REAL, DIMENSION(:), ALLOCATABLE :: ZYHAT
+TYPE(TFIELDDATA) :: TZFIELD
!-------------------------------------------------------------------------------
!
!* 1. Read orography in the file
@@ -123,7 +141,18 @@ ALLOCATE(ZSMOOTH_ZS(IIU,IJU))
ALLOCATE(ZFINE_ZS(IIU,IJU))
ALLOCATE(ZSLEVE_ZS(IIU,IJU))
ALLOCATE(ZMASK(IIU,IJU))
-!
+ALLOCATE(ZSLOPEX(IIU,IJU))
+ALLOCATE(ZSLOPEY(IIU,IJU))
+ALLOCATE(ZSLOPEMX(IIU,IJU))
+ALLOCATE(ZSLOPEMY(IIU,IJU))
+ALLOCATE(ZCOFIL(IIU,IJU))
+ALLOCATE(ZSMOOTH_ZSINI(IIU,IJU))
+ALLOCATE(ZXHAT(IIU))
+ALLOCATE(ZYHAT(IJU))
+!
+CALL IO_READ_FIELD(TPFILE,'XHAT',ZXHAT)
+CALL IO_READ_FIELD(TPFILE,'YHAT',ZYHAT)
+
!PW: bug/TODO: read a field in a file opened in WRITE mode
!There is a test in IO_READ_FIELD_BYFIELD_X2 to allow this if TPFILE%CMODE='LFICDF4'
CALL IO_READ_FIELD(TPFILE,'ZS',ZZS)
@@ -162,6 +191,8 @@ ZSMOOTH_ZS = ZZS
!
CALL ADD2DFIELD_ll(THALO_ll,ZZS)
CALL ADD2DFIELD_ll(THALO_ll,ZSMOOTH_ZS)
+CALL ADD2DFIELD_ll(THALO_ll,ZSLOPEX)
+CALL ADD2DFIELD_ll(THALO_ll,ZSLOPEY)
!
CALL UPDATE_HALO_ll(THALO_ll,INFO_ll)
@@ -192,7 +223,6 @@ CALL UPDATE_HALO_ll(THALO_ll,INFO_ll)
IF (JN==KSLEVE) ZSLEVE_ZS= ZSMOOTH_ZS
ENDDO
!
-CALL CLEANLIST_ll(THALO_ll)
!-------------------------------------------------------------------------------
!
!* 3. Case of uniform smooth orography prescribed
@@ -205,7 +235,142 @@ IF (PSMOOTH_ZS /= XUNDEF) THEN
END IF
!-------------------------------------------------------------------------------
!
-!* 4. writes smoothed orographies in the file
+!* 4. Filtering very high slopes for steep orography
+! -------------------------------------------
+!
+IF(OHSLOP) THEN
+ ZSMOOTH_ZSINI=ZFINE_ZS
+ ZZS=ZFINE_ZS
+ ZSLOPEMX=0.
+ ZSLOPEMY=0.
+ ZSLOPEX=0.
+ ZSLOPEY=0.
+ !
+ DO JN=1, KLOCZSFILTER
+ ZCOFIL=0.
+ ZSLOPEX=0.
+ ZSLOPEY=0.
+ ZSLOPEMX=0.
+ ZSLOPEMY=0.
+ !
+ !Slope calculation along Y at flux and mass point
+ DO JI=1,IIU
+ DO JJ=2,IJU-1
+ ZSLOPEY(JI,JJ) = (ZZS(JI,JJ)-ZZS(JI,JJ-1))/(0.5*(ZYHAT(JJ+1)-ZYHAT(JJ-1)))
+ END DO
+ END DO
+ !
+ CALL UPDATE_HALO_ll(THALO_ll,INFO_ll)
+ !
+ DO JI=1,IIU
+ DO JJ=1,IJU-1
+ ZSLOPEMY(JI,JJ) = MAX(ABS(ZSLOPEY(JI,JJ)),ABS(ZSLOPEY(JI,JJ+1)))
+ END DO
+ END DO
+ !
+ !Slope calculation along X at flux and mass point
+ DO JJ=1,IJU
+ DO JI=2,IIU-1
+ ZSLOPEX(JI,JJ) = (ZZS(JI,JJ)-ZZS(JI-1,JJ))/(0.5*(ZXHAT(JI+1)-ZXHAT(JI-1)))
+ END DO
+ END DO
+ !
+ CALL UPDATE_HALO_ll(THALO_ll,INFO_ll)
+ !
+ DO JJ=1,IJU
+ DO JI=1,IIU-1
+ ZSLOPEMX(JI,JJ) = MAX(ABS(ZSLOPEX(JI+1,JJ)),ABS(ZSLOPEX(JI,JJ)))
+ END DO
+ END DO
+ !
+ !Filtering coefficient
+ DO JI=1,IIU
+ DO JJ=1,IJU
+ IF(ZSLOPEMX(JI,JJ)>=PHSLOP .OR. ZSLOPEMY(JI,JJ)>=PHSLOP) THEN
+ ZCOFIL(JI,JJ) = 0.25 !Fixed coefficient
+ END IF
+ END DO
+ END DO
+ !
+ !Filtering
+ DO JJ = 1,IJU
+ DO JI = 1,IIU
+ JIP = MIN(JI+1,IIU)
+ JIM = MAX(JI-1,1 )
+ JJP = MIN(JJ+1,IJU)
+ JJM = MAX(JJ-1,1 )
+ ZSMOOTH_ZS(JI,JJ) = ZZS(JI,JJ) &
+ + ZCOFIL(JI,JJ)* ZMASK(JI,JJ) &
+ * ( ZMASK(JIM,JJ) * ZZS(JIM,JJ) &
+ + ZMASK(JIP,JJ) * ZZS(JIP,JJ) &
+ + ZMASK(JI,JJM) * ZZS(JI,JJM) &
+ + ZMASK(JI,JJP) * ZZS(JI,JJP) &
+ - ( ZMASK(JIM,JJ) &
+ +ZMASK(JIP,JJ) &
+ +ZMASK(JI,JJM) &
+ +ZMASK(JI,JJP) ) * ZZS(JI,JJ) )
+ ENDDO
+ ENDDO
+ !
+ CALL UPDATE_HALO_ll(THALO_ll,INFO_ll)
+ ZZS=ZSMOOTH_ZS
+ END DO
+ !
+ ZFINE_ZS=ZSMOOTH_ZS
+ !
+ !Filtered slopes computed again for output
+ DO JI=1,IIU
+ DO JJ=2,IJU-1
+ ZSLOPEY(JI,JJ) = (ZZS(JI,JJ)-ZZS(JI,JJ-1))/(0.5*(ZYHAT(JJ+1)-ZYHAT(JJ-1)))
+ END DO
+ END DO
+ !
+ DO JJ=1,IJU
+ DO JI=2,IIU-1
+ ZSLOPEX(JI,JJ) = (ZZS(JI,JJ)-ZZS(JI-1,JJ))/(0.5*(ZXHAT(JI+1)-ZXHAT(JI-1)))
+ END DO
+ END DO
+ !
+ ! Writes filtred orography and slopes along i and j
+ TZFIELD%CMNHNAME = 'ZSLOPEX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'ZSLOPEX'
+ TZFIELD%CUNITS = ''
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'orography slope along x'
+ TZFIELD%NGRID = 4
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 2
+ TZFIELD%LTIMEDEP = .FALSE.
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZSLOPEX)
+ !
+ TZFIELD%CMNHNAME = 'ZSLOPEY'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'ZSLOPEY'
+ TZFIELD%CUNITS = ''
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'orography slope along y'
+ TZFIELD%NGRID = 4
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 2
+ TZFIELD%LTIMEDEP = .FALSE.
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZSLOPEY)
+ !
+ TZFIELD%CMNHNAME = 'ZS_FILTR'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'ZS_FILTR'
+ TZFIELD%CUNITS = 'm'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'filtred orography'
+ TZFIELD%NGRID = 4
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 2
+ TZFIELD%LTIMEDEP = .FALSE.
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZSMOOTH_ZSINI-ZFINE_ZS)
+END IF
+!-------------------------------------------------------------------------------
+!
+!* 5. writes smoothed orographies in the file
! ---------------------------------------
!
!
@@ -217,6 +382,14 @@ DEALLOCATE(ZFINE_ZS)
DEALLOCATE(ZSMOOTH_ZS)
DEALLOCATE(ZSLEVE_ZS)
DEALLOCATE(ZMASK)
+DEALLOCATE(ZSLOPEX)
+DEALLOCATE(ZSLOPEY)
+DEALLOCATE(ZSLOPEMX)
+DEALLOCATE(ZSLOPEMY)
+DEALLOCATE(ZCOFIL)
+DEALLOCATE(ZSMOOTH_ZSINI)
+!
+CALL CLEANLIST_ll(THALO_ll)
!
!-------------------------------------------------------------------------------
!
diff --git a/src/SURFEX/init_surf_atmn.F90 b/src/SURFEX/init_surf_atmn.F90
index 6e454467d9b2873f4ff877cc14a5906572cfb79f..2c6c776dc8c4efdcea39793fc739980dada1926e 100644
--- a/src/SURFEX/init_surf_atmn.F90
+++ b/src/SURFEX/init_surf_atmn.F90
@@ -55,6 +55,7 @@ SUBROUTINE INIT_SURF_ATM_n (YSC, HPROGRAM,HINIT, OLAND_USE, &
!! R. Séférian 03/2014 Adding decoupling between CO2 seen by photosynthesis and radiative CO2
!! M.Leriche & V. Masson 05/16 bug in write emis fields for nest
!! (P.Tulet & M.Leriche) 06/2016 add MEGAN coupling
+!! J.Escoabr 01/2019 integrate bypass fo albedo pb > 1.0 from Florian Pantillon (Sep 2011)
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -219,6 +220,8 @@ REAL :: XTIME0
INTEGER :: ISIZE_FULL
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+INTEGER :: JJ
!-------------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('INIT_SURF_ATM_N',0,ZHOOK_HANDLE)
@@ -663,6 +666,19 @@ IF (SIZE(PTSURF)>0) &
!
DEALLOCATE(ZFRAC_TILE)
!
+! MODIF FP SEP 2011
+DO JJ=1,KI
+ IF (PDIR_ALB(JJ,1)>1.) THEN
+ WRITE (*,*) 'JJ', JJ
+ WRITE (*,*) 'PDIR_ALB', PDIR_ALB(JJ,:)
+ WRITE (*,*) 'PSCA_ALB', PSCA_ALB(JJ,:)
+ WRITE (*,*) 'PEMIS', PEMIS(JJ)
+ WRITE (*,*) 'PTSRAD', PTSRAD(JJ)
+ PDIR_ALB(JJ,:) = 0.5
+ PSCA_ALB(JJ,:) = 0.5
+ END IF
+END DO
+! END MODIF FP SEP 2011
!-------------------------------------------------------------------------------
!==============================================================================
IF (LHOOK) CALL DR_HOOK('INIT_SURF_ATM_N',1,ZHOOK_HANDLE)